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Abstract. We use the Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 9 
(DR9) to detect and measure the position of the Baryonic Acoustic Oscillation (BAO) fea- 
ture in the three-dimensional correlation function in the Lyman-o: forest flux fluctuations 
at a redshift Zeff = 2.4. The feature is clearly detected at signiflcance between 3 and 5 
sigma (depending on the broadband model and method of error covariance matrix estima- 
tion) and is consistent with predictions of the standard ACDM model. We assess the biases 
in our method, stability of the error covariance matrix and possible systematic effects. We 
flt the resulting correlation function with several models that decouple the broadband and 
acoustic scale information. For an isotropic dilation factor, we measure 100 x (a[so — 1) = 
— 1.6^2 -4 1 -6 8 (stat.) zbl.O (syst.) (multiple statistical errors denote 1,2 and 3 sigma confl- 
dence limits) with respect to the acoustic scale in the flducial cosmological model (flat ACDM 
with — 0.27, h — 0.7). When fltting separately for the radial and transversal dilation fac- 
tors we flnd marginalised constraints 100 x {a\\ — 1) = —1.3^3 3 7 ^lo 2 (stat.) ±2.0 (syst.) 
and 100 x [a^ — 1) = — 2.2^71 (stat.) ±3.0 (syst.). The dilation factor measurements are 
signiflcantly correlated with cross-correlation coeflicient of ^ —0.55. Errors become signifl- 
cantly non-Gaussian for deviations over 3 standard deviations from best flt value. Because 
of the data cuts and analysis method, these measurements give tighter constraints than a 
previous BAO analysis of the BOSS DR9 Lyman-a forest sample, providing an important 
consistency test of the standard cosmological model in a new redshift regime. 

Keywords: cosmology, Lyman-o^ forest, large scale structure, dark energy 
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1 Introduction 



Discovery of the accelerated expansion of the Universe at the end of 20* century was one 
of the largest surprises in cosmology [1, 2]. Since then, much effort has gone into measuring 
the properties of the dark energy that is believed to drive the acceleration of the Universe 
using a variety of techniques (see extensive review [3] and references therein). Measurements 
that directly probe the expansion history of the Universe reveal how the energy density of its 
constitutive components changes with expansion and thus allow inferences to be made about 
their nature. The two important techniques in this field are the type la supernovae (SN la) 
and baryonic acoustic oscillations (BAO) [4-9]. The former rely on the fact that SNIa are 
standarizable candles. Measuring supernovae fluxes and redshifts thus provide information 
about the luminosity distance as a function of redshift. The BAO technique, which is the 
focus of this paper, relies on the fact that the baryonic features in the correlation properties 
of tracers of the large scale structure of the Universe can act as a comoving standard ruler. 

BAO are acoustic oscillations in the primordial plasma before decoupling of the bary- 
onic matter and radiation. These oscillations imprinted a characteristic scale into the corre- 
lation properties of dark matter, which exhibit themselves as a distinct oscillatory pattern 
in the power spectrum of fluctuations or equivalently a characteristic peak in the correlation 
function. The scales involved are large 150 Mpc) so the feature remains in the weakly 
non-linear regime even today. Any tracer of the large-scale structure can measure this feature 
and use it as a standard ruler to infer the distance to the tracer in question (when the ruler 
is used transversely) or the expansion rate (when the ruler is used radially). 

The main attraction of BAO technique is that one is measuring the position of a peak 
and that many systematic effects are blind about the preferred scale and thus influence 
the two-point function only in broadband sense. Such systematic effects can be efficiently 
marginalized out with little or no signal loss. The caveat in the case of the Lyman-a forest 
analysis is that systematic effects that appear at certain pairs of wavelengths can produce a 
peak-like feature in the correlation function in the radial direction. 

Traditionally, the tracer used to measure BAO consisted of galaxies, and measuring 
BAO with galaxies is a mature method [8, 9]. However, the current measurements of the 
BAO with galaxies are restricted to redshifts z < 1, due to difficulty in measuring redshifts 
of sufficient numbers of objects efficiently. Even the proposed future BAO projects such as 
BigBOSS [10] will measure the BAO only to redshifts of less than z ^ 1.5 — 2 using galaxies as 
a tracer of cosmic structure. However, there is strong scientific motivation for measuring the 
expansion history of the Universe at higher redshifts. Dark energy suffers from many tuning 
problems that can be alleviated by introducing an early dark energy component [11-16]. 
Such a component can only be measured by a technique sensitive to the expansion history 
at high redshift. In this paper we present such a measurement using the Lyman-a forest as 
a tracer of dark-matter fluctuations at z ^ 2.4. 

This "Lyman-a forest " denotes the absorption features in the spectra of distant quasars 
blue- ward of the Lyman-o; emission line [17]. These features arise because the light from 
a quasar is resonantly scattered by the presence of neutral hydrogen in the intergalactic 
medium. Since the quasar light is constantly red-shifting, hydrogen at different redshifts 
absorb at different observed wavelength in the quasar spectrum. The amount of scattering 
reflects the local density of neutral hydrogen. Neutral hydrogen is believed to be in photo- 
ionization equilibrium and therefore the transmitted flux fraction at position x is given by 
[18] 
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F(x) = exp[-r(x)] ~ Cexp [-^(1 + Sb{^)r] , (1.1) 

where r is optical depth for Lyman-o; absorption at position x. This optical depth scales 
roughly as a power law in the baryon over-density 6^,. The coefficient p would be two for 
an idealized two-body process, but is in practice closer to ^ 1.8 due to the temperature 
dependence of the recombination rate. The Lyman-o: forest thus measures some highly 
non-linear but at the same time very local transformation of the density field. Any such 
tracer will, on sufficiently large scales, trace the underlying dark matter field in the observed 
redshift-space as 

SF{k) = ib + b,f^i^)Sm{k), (1.2) 

where ^^(k) are the flux and matter over-densities in Fourier space respectively and 6jn is 
the matter over-density in the Fourier space. The quantity /i = fc||/|k| is the cosine of the 
angle of a k vector with respect to the line of sight and the parameter / = ding /din a is the 
logarithmic growth factor. In the peak-background split picture, the two bias parameters 
can be derived as a response of the smoothed transmitted flux fraction field Fg to a large 
scale change in either overdensity 5/ or peculiar velocity gradient r]i = —H~'^dv^^/dr^^ [19]: 

= (1.4) 

dm 

Both bias parameters (6 and by) can in principle, be determined from numerical simu- 
lations of the intergalactic medium, but they can also be determined from the data as free 
parameters to be fit. The dimensionless redshift-space distortion parameter /3 = fK/b deter- 
mines the strength of distortions [20]. Conservation of the number of galaxies requires that 
by = 1 in that case and hence /3 is completely specified by the value of / and density bias 
b. Unfortunately, there is no such conservation law for transmitted flux fraction and hence 
the two parameters are truly independent for the Lyman-a forest 2-point function. In the 
rest of this paper we will work with density bias parameter b and redshift-space distortion 
parameter /3 and treat them as free parameters to be constrained by the data. 

Numerical simulations have shown that for plausible scenarios, the Lyman-a forest is 
well in the linear biasing regime of equation (1.2) at scales relevant for BAO. Fluctuations 
in the radiation field can, if sufficiently large, invalidate the equation (1.2), but it has been 
shown that these variations are unlikely to produce a sharp feature that could be mistaken 
for the acoustic feature [21, 22]. Therefore, measuring three-dimensional (3D) correlations in 
the flux fluctuations of the Lyman-o; forest provides an accurate method for detecting BAO 
correlations [23-25]. 

Using the Lyman-a forest to measure the 3D structure of the Universe became possible 
with the advent of the Baryonic Oscillation Spectroscopic Survey (BOSS, [26, 27]), a part of 
the Sloan Digital Sky Survey III (SDSS-III) experiment [28-33]. This survey was the first to 
have a sufficiently high density of quasars to measure correlations on truly cosmological scales 
[34] . This paper was also ffist to point out difficulties in measuring the large-scale fluctuations 
in the Lyman-a forest data; and many design decisions in this paper were chosen to alleviate 
them. The BAO in the Lyman-a forest was first detected in [35] using the same dataset as 
this work. They found a dilation scale aiso = 1.01 ± 0.03. Both works were done within 
the same working group in the BOSS collaboration, but use different and an independently 
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developed pipelines. We assess these differences more thoroughly in Section 5.4. We refer 
the reader to these publications for a more thorough review of the Lyman-a forest physics 
and the BOSS experiment. 

This paper is focused on the methodology and results of measuring the 3D correlation 
function from the flux transmission fluctuations. The companion paper [36] focuses on the 
methodology of constraining the BAO parameters. The final results on the BAO parameters 
are presented in this paper. 

This paper is dedicated purely to detecting the BAO feature in the correlation function of 
Lyman-a forest fluctuations as measured by BOSS Data Release 9 [37]. We do not attempt 
to provide robust measurements of bias and /3 parameters; we therefore adopt aggressive 
inclusion of data that is tuned to maximize the likelihood of a significant BAO peak detection. 
We stress that in the early stages of analysis, the peak between 100-130/i~^Mpc was blinded 
and these cuts were decided on and fixed before "opening the box" . No cuts were changed 
after that, but the data analysis code has been improved. 

The paper is structured as follows. In section 2 we describe our data samples and 
basic selection criteria. In section 3 we explain the complicated process of measuring the 
3D correlation function using an optimal estimator, its errors and biases, while presenting 
intermediate results. Section 5 discusses application of the method to the real and synthetic 
data and discusses potential errors. Finally, we convert our many noisy measurements into 
a plot of the BAO bump in Section 6. We conclude in Section 7. The mathematical details 
of optimal estimators are elaborated in Appendices. 

2 Instrument, Data & Synthetic data 
2.1 DR9Q Sample 

In this work we use BOSS quasars from the Data Release 9 (DR9, [37]) sample in order to 
facilitate reproducibility by outside investigators. The quasar target selection for the DR9 
sample of BOSS observations is described in detail in [38] and [39]. A similar dataset is 
released in [40] including continua for normalization and default recommendations for data 
masking. 

The main approach in this work has been to aggressively use as much data as possible, 
to maximize the likelihood of a significant detection of the BAO peak and the accuracy of the 
determination of its position. The argument that we are making, and which underpins the 
robustness of many BAO measurements, is that systematic effects that affect the measured 
two-point statistics will in general produce a broadband contribution and thus not affect 
the position of the peak itself, since broadband contributions are marginalized out in the 
process of fitting the BAO position. At the same time, however, adding poorly understood 
data might lead to an increase in the noise that might outweigh the benefit of additional 
information. We settled on a set of cuts described below. 

We select quasars from the DR9Q version of the Value Added Catalog (VAC) of the 
French Participation Group (FPG) visual inspections of the BOSS spectra [41]. The position 
of the quasars on the sky is plotted in figure 1. The total number of quasars in our analysis 
after the cuts performed below is 58,227. The covered area is optimized through survey 
strategy considerations and will become more uniform as the survey approaches completion 
in 2014. 

We limit our analysis to quasars with redshifts 2.1 < z < 3.5. Throughout this work 
we use the visual redshifts (Z_VI) from the DR9Q. The low number density of quasars at 
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Figure 1. Distribution of 58,227 quasars used in this work on the sky in J2000 equatorial coordinates, 
shifted by 180 deg to more clearly show the main survey area. 

redshifts larger than 3.5 renders them not useful for measuring the BAO signal, which is the 
main goal of this paper. We show the distribution of quasar redshifts and the Lyman-a forest 
pixels in Figure 2. 

For repeated observations of a given object the DR9Q catalog contains all the hsted 
observations, as well as the observation deemed to be optimal. We performed the analysis 
either by restricting the sample to the subset of best observations as well as using all ob- 
servations. In the latter case, we simply co-add the observations on the nearest pixel point 
using the pipeline-reported inverse variances. We do not use the noise correction (discussed 
later) when performing this co-add, but since the correction is a function of wavelength it 
would rescale all contributions equally and thus have no effect when stacking in observed 
wavelength. 

We apply several cuts to the data. First, we treat quasars that are considered Broad 
Absorption Line (BAL) quasars differently. The DR9Q catalog contains a measured value 
of balnicity index [42]. We remove quasars with BI_CIV > 2000 km/s. We additionally drop 
quasars that DR9Q visual inspection flagged as highly unusual (BI_CIV = —1). 

Second, we mask regions around bright sky lines. The mask is obtained by the following 
procedure: We subtract the sky model from the sky fibers. The model's prediction at the 
position of the sky fiber does not match the spectrum of the sky fiber, since the model fits a 
smooth polynomial model over all sky fibers and thus sky subtracted sky fibers do not have 
zero flux. We stack all sky-subtracted sky fibers and measure the root mean square defined 
over ±25 standard BOSS pixels (A log^^Q A = 10~^) boxcar of this stack. Next we mask all 
pixels that exceed this r.m.s. by a factor of 1.25. The masked pixels are excluded from the 
r.m.s. measurement and the process is iterated until the mask converges. This mask excludes 
very few pixels in the blue end of spectrograph, simply because there are not many sky-line 
in that wavelength region. 
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Figure 2. Distribution of quasar redshifts (red, solid) and Lyman-a forest pixels (blue, dotted). 

Third, we exclude portions of the data affected by the Damped Lyman-a (DLA) sys- 
tems using the so-called concordance DLA sample [43]. Briefly, this dataset was established 
by three groups within the BOSS Lyman-a forest working group, two of which developed 
algorithms for detecting DLAs in data automatically and the third inspected all quasars 
visually. All three groups have provided catalogs of systems. One of these three samples is 
published in [44]. The concordance catalog DLAs were identified by at least two of the three 
groups. We do not use data within 1.5 equivalent widths from the center of a DLA. 

When measuring the correlation in the Lyman-a forest, we use a wide region of the 
forest, spanning 1036-1210 A for most quasars and 1036-1085 A for the remaining BAL 
quasars (0 < Bl < 2000). We also exclude pixels with observed wavelength less than 3600 A, 
beyond which the signal-to-noise of the spectrograph and spectro-photometric calibration 
become very poor. 

We assume that the measurement noise does not correlate between pixels, but we apply 
a correction to the BOSS spectral extraction pipeline estimates of the noise. To derive this 
correction, we fit a cubic polynomial to the smooth region of quasar spectra between rest- 
frame 1420 and 1510 A and compared the fit residuals with the pipeline error estimates. We 
find that the co-added spectra have an observer-frame wavelength dependent mis-estimation 
of the noise, which follows the square root of the ratio of pixel sizes between the co-added 
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and individual exposure spectra. This value is a — 10% under-estimate of the noise variance 
below observer-frame ^4750 A and a — 15% over-estimate above ^4750 A; this correction 
is not valid beyond ^6000 A where both arms of the BOSS spectrograph contribute to the 
spectrum. We correct for this mis-calibration by adjusting the pipeline noise by the square 
root of the co-add-to-individual exposure pixel size ratio. 

2.2 Synthetic data 

Data analysis procedures in this work were tested on the synthetic data that was generated 
in nearly the same manner as for the previous work [34] and we briefly review it here for 
completeness. 

We want to generate a transmitted flux- field with desired two-point function and approx- 
imately correct probability distribution function. For a general transformation F = T{6g)^ 
where Sq is a Gaussian auxiliary fleld, the 2-point correlation functions of the flux fleld, 

= (5i?(x)5i?(x + v)), and that of the auxiliary Gaussian fleld, ^(|v|) = (5^(x)5^(x + v)), 
can be related by 



= / dFi [ dF2p{Fi,F2)FiF2 
Jo Jo 



d6i / dS2p{Si,S2)T{Si)T{S2) 



OO POO 2(l-e2(ri2)) 

dS^ / d52—^=^^m)T{S2) . (2.1) 



■OO 



We assume that r = log F is log- normally distributed so that specifying the mean trans- 
mitted flux fraction F, in addition to the desired target power spectrum completely specifles 
the required properties of the auxiliary Gaussian fleld. This inversion is done numerically. 
The Gaussian fleld is initially generated assuming parallel Lyman-a forest lines of sight and 
stationary (non-evolving) fleld: the effects of non-parallel lines of sight and fleld evolution 
are then obtained by interpolating between redshifts. The assumption of parallel lines of 
sight and stationary flelds signiflcantly simplifles the problem. Decomposing the line of sight 
flux fluctuations of a quasar q into its Fourier components 5/c(^, ^id), the cross-correlation 
vanishes for different Fourier modes: 

{Sk{q, hMq', k[o)) = S{k - k')P:,{hD,r±{q, q')). (2.2) 

This allows one to generate the fleld one parallel Fourier mode at a time. For each such- 
mode, one thus needs to generate a correlated Gaussian fleld vector of the size Nq, where Nq 
is the number of quasars. We perform this using standard methods for generating correlated 
Gaussian flelds with a desired covariance matrix by Cholesky decomposing the matrix and 
multiplying with a vector of uncorrelated unit- variance Gaussian fleld. Technical aspects of 
synthethic data generation and tests that demonstrate the validity of the generated flelds are 
described in great detail in [45]. 

Our present implementation of synthethic data generation is still numerically extremely 
demanding. To be able to generate a full synethetic dataset, we split it into four geometrically 
self-contained "chunks", ignoring correlations across chunk borders. We have tested explicitly 
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that this approximation introduces no bias. We have thus generated fifteeen synthethic 
reahzations of our dataset. These synethetic datasets were also used in [35]. 

Compared to [34] , mock data were created using shghtly different bias parameters b = 
—0.14 and /3 = 1.4. There were made more reahstic in this work by adding several effects 
that are present in the spectrograph and were not accounted for in the previous work. We 
add noise to the mocks based upon an object-by-object fit to the observed photon noise in 
the real data. We then generate a noisy and mis-calibrated estimate of the true noise to 
model the pipeline reported noise and use that noise to generate a noise realization for the 
analysis. The mock spectra fluxes were multiplied by a polynomial in A to match the real 
data to simulate QSO flux mis-calibration and a small amount of sky spectrum is added to 
model a sky subtraction bias in the BOSS data. These additions will be described in greater 
detail in a future publication [46]. In one respect, the synthetic data were less realistic that 
the similar synthetic datasets used in [34], namely, we did not contaminate the data with 
associated forest metals in this work. 

3 Data Analysis 

In this section we describe the steps taken in the data analysis. The same procedure is 
adopted both to the real and the synthetic data described in section 2.2. 
The overview of data analysis is as follows: 

• Continuum fitting. (Section 3.1) First we perform the basic continuum fitting, where 
the pixels are treated as independent. Continuum fitting can be described as a process 
of fitting a model for what the spectrum of a given quasar would look like in a com- 
pletely homogeneous universe, i.e. there would be no forest, except the mean hydrogen 
absorption. We fit for the mean continuum, the redshift dependence of the mean trans- 
mission in the forest and the intrinsic variance including its redshift dependence. We 
also fit for the per-quasar amplitude and slope across the forest. 

• One- dimensional Pixel Power Spectrum measurement. (Section 3.2) In the second step 
we generalize the variance in the pixels to the full one-dimensional (ID) power spectrum 
of flux fluctuations as a function of redshift. We refit the individual quasar forest 
amplitudes and effective spectral slopes, the mean transmission as a function of redshift 
and the effective ID power spectrum. This power spectrum is not appropriate for 
cosmological analysis, as it does not deconvolve the instrument PSF and the instrument 
noise is treated coarsely. This co-variance is only required for the weighting in the 
measurement of the 3D correlation function. 

• Estimation of the three-dimensional correlation function and its covariance matrix. 
(Section 3.3) The next step is estimation of the 3D correlation function. This is done 
by weighting quasar data with a per-quasar inverse covariance matrix and performing 
optimal quadratic estimation of the correlation function. This analysis is done in a 
coordinate system that is composed exclusively of observable quantities, i.e. difference 
in logarithm of the observed wavelength, angular separations on sky and the pixel pair 
redshift (which is just a proxy for the mean redshift of a given pixel pair). We use 
several methods to estimate the covariance matrix coming from the optimal estimator. 

• Estimation of BAO position and significance. (Section 4) Finally, we measure the 
position of the BAO peak from the measured 3D correlation function. The uncertainties 
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on the BAO parameters of interest using different covariance matrices and broadband 
models. 



3.1 Continuum fitting 

A reasonably general model of the observed flux can be written as follows 
fiq, \o) = PSF(g, Xo) * [A{q, Xr)Fi\o)CiXr)T{Xo){l + Spiq, A^)) + 5(Ao)] + e(Q, A^). (3.1) 

Here we use the index q to denote a particular quasar (probing a direction on the sky), 
Xo is the observed wavelength and = Xq/ + is the rest-frame wavelength of the quasar. 
Since hydrogen absorbing through Lyman-o: at redshift z produces a decrement at observed 
Ao = (1 + z)Xa^ we can use A^ interchangeably with absorption redshift. The PSF(g, A^)^ 
term denotes convolution, i.e., smoothing by the point-spread function of the spectrograph 
and the final term in equation (3.1) is the spectrograph noise contribution. The pipeline 
provides an estimate of the pixel noise under the assumption that it is independent from 
pixel to pixel (but we correct for this as described earlier). 

The mean unabsorbed quasar spectrum, also referred to as the continuum, is given 
by C{Xr) and we absorb the diversity in both the quasar amplitude and quasar shape into 
A(q^ Xr). In this work, we model A{q^ A^) as linear function of log A 

A{q, Xr) = A{q, Xoi) + (A(g, A,,) - A{q, Xa)) ^""^ ~ , (3.2) 

log Xoh - log Xoi 

where XqH and Xoi are the highest and lowest wavelengths in the forest. 

The F{Xo) term in Equation 3.1 is the mean quasar transmission in the forest and can 
be distinguished from the residual variations in the spectrograph transmissivity, T{Xo), by 
the virtue that the former is unity outside the forest. Finally, the ^(A^, A^) term describes an 
additive systematic component, likely to be dominated by the sky. The 6f term inside brack- 
ets are the fluctuations in the transmission, whose correlation properties we are attempting 
to infer from the data. 

With expediency in mind, we make several simplifying assumptions. First, we ignore 
the point-spread function of the spectrograph, since it affects exclusively small scales and 
can be completely neglected for purely 3D analysis. Second, we set the sky contribution 
^(A^, Xo) = and the spectrograph transmissivity T{Xo) = 1. This choice means that C(A^) 
and F{Xo) are really just parametrization of continuum model under the assumption that it 
is factorisable as C{Xr)F{Xo) after the individual quasar mean luminosity and spectral slope 
have been fitted (Equation (3.2)). The fitted F{Xo) will also absorb the contribution from 
the variation of spectrograph properties as a function of wavelength which have not properly 
been calibrated by the pipeline. 

Finally, in the initial parts of the fitting process, we model fluctuations as Gaussian and 
uncorrelated. Therefore, our likelihood model in the first step of continuum fitting can be 
written as 



logL = 5] 



{fiq, Xo) - A{q, Xo)C{Xr)FiXo)f log[iV(g, A^) + a(A„ 



2iN{q,Xo)+a^{Xo)) 



(3.3) 



where the sum is over all forest pixels in all quasars, denotes the intrinsic variance of 
the flux transmissions (i.e., variance in the cosmological 6^ field) and N = (e^) the pipeline 
noise. 
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We measure two amplitude parameters for each quasar {A(q, X^i) and A(q^ Xoh)^ defining 
amplitude and slope), C{Xr) in 95 bins between = 1036A and = 1210A and further 
167 bins extending to A^ = 1600 A. The intrinsic variance and F are measured in 16 bins 
between redshifts 1.95 and 3.45 (linearly interpolating between bins for the case of F). 

This approach contains a considerable number of parameters. We maximize the likeli- 
hood by brute-force maximization, employing a Newton-Raphson algorithm with numerical 
derivatives. We do not attempt to derive errors on these parameters as these errors will not 
dominate the errors on the final quantities that we are attempting to derive. 

To gracefully approach a good solution, we start by fitting just C, F, intrinsic variances 
and constant amplitudes across rest wavelength 1036 - I6OOA. Next we focus on the forest 
region and refit individual amplitudes and slopes together with refitting variances. As a 
third step we refit variances and slopes together with measuring the ID power spectrum as 
described in the next section. 

Figures 3 and 4 show the mean continuum and F for data and synthetic data. The 
measurement of F shows the expected structure: there is more absorption at higher redshifts. 
However, we also see structure at the position of rest-frame Balmer lines. This structure is 
an artifact of the pipeline and arises due to imperfect interpolation of masked Balmer regions 
in the calibration stars. These features are essentially completely absorbed into F and we 
do not see any residual structure associated with Balmer lines in other tests. We also see 
features at the positions of CaH and CaK absorptions. 

In Figure 5 we show the distribution of amplitudes and effective spectral slopes as fitted 
to the data and to the synthetic data. The histogram of amplitudes for the data matches 
that of the synthetic data very well. This result arises essentially by construction, since 
our synthetic quasars are matched in flux to real quasars. The right-hand panel considers 
the spectral slope. The slope is defined so that it is ±1 if any one of the two individual 
fitted amplitudes at two wavelengths is zero. We fit for the slope at two fixed observed 
wavelengths at the edge of the forest even though there might be no data at those edges. 
Therefore the amplitude measurement can be quite uncertain and produce a value outside a 
physically acceptable range. In general, the diversity of spectra in synthetic data is a good 
representation of that in real data, although, as expected the real data are slightly more 
diverse. 

Our model for the continuum in Equation 3.1 is not perfect and there will be systematic 
residuals between true 6f and our estimates in addition to noise. A useful check for these 
residuals is to stack our estimates of 6f in bins of observed wavelength (or redshift): for 
unbiased estimates these should vanish, but in practice this is not exactly the case. In Figure 6 
we show the residuals from our continuum fitting which are obtained by an inverse covariance 
weighted average of all data in bins of observed wavelength. These residuals are removed 
from the data before measuring the 3D correlation function, but they essentially imply a 
calibration uncertainty on the final correlation function of a couple of percent (although 
this is irrelevant for BAO fits). Since we are fitting for the amplitude of quasars in flux 
and then dividing by a noisy estimate to derive we expect that the mean of 6f will 
not be zero (essentially because 1/ (A) ^ {1/A)). We also notice that these residuals are 
considerably worse for the synthetic data. This result arises partially because the synthetic 
data is somewhat noisier, especially at higher redshift (see the next section) and partially 
because addition of some systematic effects in the synthetic data might be introducing more 
contamination than exists in the real data. Most importantly, however, we see no effects of 
structure at the position of Balmer lines, indicating that these features are nearly perfectly 
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Figure 3. Fitted mean continuum, C(Arest) in arbitrary units per unit wavelength. Solid red curve 
corresponds to data, while the dashed green line is for synthetic data. 



removed with the F fit. The caveat is, of course, that one could have structure associated 
with Balmer features that would be, for example, correlated from plate to plate and produce 
contaminating signal while still averaging to zero. This concern is discussed further in Section 
5.2. 

3.2 One-dimensional Pixel Power Spectrum measurement 

As the final step in the continuum fitting process, we measure the ID power spectrum of 
pixels. This power spectrum enables a nearly optimal weighting to be used when measuring 
the 3D two-point function. This step essentially generalizes the model of the covariance of 
the flux transmission A^) from uncorrelated pixels to pixels correlated across a single 
quasar. 

(%,Ao)%,A:,)> =e(^eff,AlogAo) = PiD(Zeff,fc)e^'='^l°sAo^^ (3.4) 

In practice, we measure the power spectrum in flat power bands in k (i.e., the ID power 
spectrum is approximated as a set of top-hat bins) and interpolated in redshift. For any pair 
of pixels, the effective redshift Zgff is taken to be given by the geometrical mean of (1 + z) 
of two pixels (or equivalently geometrical mean of their Arest mean of their log Arest)- The 
distance is measured in Alog A^ = | log A^ — log A^l converted to the standard units of km/s 
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Figure 4. Fitted mean transmission in the forest {F) with arbitrary normahzation. Sohd red curve 
corresponds to data, while the dashed green hne is for synthetic data. Dotted hnes show the Balmer 
series. The blue dashed lines are positions of Galactic calcium extinction. 

(by multiplying by speed of light). The integral of flat power band can be trivially calculated 
analytically. 

To estimate the ID power spectrum we use the standard optimal quadratic estimator 
[47-49] and has been applied to Lyman-o; forest data in before in [50]. The formalism 
behind optimal quadratic estimators generalizes the signal-to-noise weighting of the data from 
scalar weights proportional to the inverse variance applied to each data-point individually, to 
multiplying the entire data- vector by the inverse of the covariance matrix. Optimal quadratic 
estimator methodology can improve the signal-to- noise of the measurement and provides good 
(but not perfect) error-estimates. We review the methodology in the Appendix A. 

We measure the ID power spectrum in 20 bins k logarithmically spaced between 10~^-^^ 
s/km (center of the first bin) and lO"-*^*^^ s/km (center of the last bin) in steps of 0.1 dex 
and 9 bins in redshift, uniformly spaced between 1.8 and 3.4. Although we have very few 
pixels with z < 2.0, they contribute to the power spectrum measurements extrapolated to 
the z = 1.8 bin (to keep the binning in redshift uniform). This results in 180 ID power 
spectrum overall bins. The lower k range for the lowest fc-bin is extended to /c = 0. 

We do not attempt to deconvolve the point spread function or carefully understand the 
effect of noise - we are only attempting to estimate the correct pixel-pixel correlations to use 
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Figure 5. The amplitude and slope distributions when fitting data and synthetic data. The left- 
hand slide plot shows the histogram of fitted amplitudes for data (red, solid) and synthetic data (green, 
dashed). The right-hand side plot shows the distribution of amplitudes and slopes for data (red points) 
and synthetic data (green points). See text for further discussion. 

in weighting in the 3D correlation estimation. For example, if noise is underestimated, then 
one measures a power spectrum that has a constant offset with respect to the true power 
spectrum, but the total signal-plus-noise power will remain unchanged. The same is true for 
the effect of the point-spread function. If we deconvolved the effect of the PSF, we would 
be putting exactly the same PSF back into the correlation matrix when making a prediction 
for the covariance matrix of a single quasar forest. In summary, in this paper, the ID power 
spectrum is treated as an intermediate product used for weighting the 3D two-point function 
estimate and is not intended to be useful as a measurement of the physical ID Lyman-a forest 
power spectrum. Any error in ID power spectrum results in a less-optimal 3D weighting, 
but will not bias the 3D measurement. 

When measuring the ID power spectrum, we iterate between measuring the ID power 
spectrum and refitting all per-quasar amplitude and slope parameters. The process converges 
in approximately five iterations. When measuring the ID power spectrum we consistently 
marginalize out modes that were fitted by the continuum fitting as described in Appendix B 
and discussed later. 

In Figure 7 we plot the ID power spectra for real and synthetic data. The data and 
the synthetic data show a good overall agreement, although they differ considerably in the 
high redshift bins. The edge redshifts are essentially extrapolations from the measured data 
and should not be taken too seriously. In addition, the data in the highest redshift bins 
come from the region close to the quasar (we use data all the way to 1210A and do not use 
quasars with z > 3.5). Finally, in this work we use only a rough correction for the pipeline 
mis-estimation of the noise, which can be thought of as a constant additive uncertainty for 
the data. The errors reported by the synthetic datasets are biased (on purpose, to model the 
pipeline errors) which we did not attempt to correct, so it is natural that the synthetic data 
will not measure just the intrinsic power spectrum. However, when doing inverse covariance 
weighting for the 3D measurement, the relevant quantity is the sum of noise and intrinsic 
power, so noise mis-estimation will result in extra (or less) intrinsic power, but will affect 
the sum only at second order. We also plot the measured values from [50], which have been 
corrected for the effect of resolution and thus do not show the beam suppression at high k 
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Figure 6. Residual estimates of 5f averaged over bin in the observed frame wavelength. Green is 
for one realization of the synthetic data, while red is for the real data. Dotted lines show the Balmer 
series that is clearly seen in the F fitting in the Figure 4. The blue dotted lines are positions of 
Galactic calcium extinction. 



values. The wiggles in the measured power spectra due to presence of Si III are clear in both 
our dataset and [50] power spectra, while absent from the synthetic data. 

3.3 Quasi-optimal estimation of the three-dimensional correlation function 

Once the ID correlation properties of quasars are measured, we proceed to measure the 3D 
correlation function. 

As in the case ID power spectrum measurement, we use the optimal quadratic estima- 
tor. Applying the estimator blindly would involve inverting matrix the size of the number of 
pixels. Since this is numerically prohibitive, we instead approximate it as a block-diagonal 
matrix, where each block corresponds to a single quasars. In other words, for the purpose 
of data-weighting, we assume quasars to be uncorrelated. This approach is optimal in the 
limit of small cross-correlations between neighbouring quasars. Technical details are out- 
lined in the Appendix A, in particular, we numerically implement equations A. 17, A. 18 and 
A. 19. This approach is always unbiased, but the resulting error-covariance matrix will ignore 
contributions to the error covariance matrix stemming from correlations between pixels in 
neighboring quasar pairs. 
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Figure 7. One-dimensional power spectra for data (solid red curve) and synthetic data (dashed green 
curve). Panels correspond to redshifts bins indicated. Noise has been subtracted, but the PSF has not 
been corrected for. Blue and black curves are measurements from [50], before and after corrections 
for background power. See text for further discussion. 

We never measure the 3D correlations using pixel-pairs that reside in the same quasar 
since continuum fitting errors are expected to be correlated within the same quasars, but 
significantly suppressed across quasars (true quasar continua are thought to be completely 
uncorrelated across quasars^). We sort pairs of quasars into SDSS plates. If a given quasar 
pair has two objects observed on two different plates, the pair is uniquely associated with 
one of the two plates. There are 817 plates in DR9 which parcels the problem of computing 

notable exception are systems with multiply lensed quasars. 
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the 3D correlation function into feasible computation chunks, but also allows a suitable base 
for bootstrapping errors. In particular, one could imagine a class of systematic errors that 
are associated with plates and such errors would be "discovered" as extra variance by the 
bootstrap algorithm. 

However, even in the approximation where quasars are decoupled, the computation of 
estimates is still too slow. Therefore, we reduce the problem by optimally reducing the 
number of pixels by a factor of three or four. This reduction is done after multiplying the 
data vector by the inverse covariance matrix as described in Appendix C. 

Our model for the 3D correlation function is parametrized in terms of completely ob- 
servable coordinates that do not assume a fiducial cosmology. In particular, each bin of the 
measured correlation function is characterized by the radial separation measured in A log A 
(where A is the observed wavelength of a given pixel) , angular separation measured in radians 
and redshift (which is just a proxy for the mean observed wavelength). In our parametrized 
model for the correlation function that we measure, we interpolate linearly in A log A and 
log(l + z). This is particularly important in the redshift direction, where bins are quite 
wide. In the angular direction, however, the bins are assumed to be uncoupled (i.e., top hat 
in shape). The reason for this choice is that it significantly simplifies the resulting Fisher 
matrix, which is zero for elements corresponding to pairs of estimates at different separations. 
This point will be elaborated upon further in Section 3.5.2. 

The main advantage of this approach is that one is truly measuring correlations without 
an assumption on cosmology. The main drawback is the large number of correlation bins 
and an inevitable smoothing of the measured correlation function, because the BAO peak 
position in (5 log A, 5) configuration space is a function of redshift. We will show later that 
this is a neghgible effect. We use three redshift bins {z = 2.0, z = 2.5 and z = 3.0), 18 
separation bins centered around 5,15,. . . 175 arc minutes and 28 bins in A log A spanning the 
relevant range with non-uniform bin spacing (0, 0.001, 0.003, 0.005 . .. 0.049, 0.059, 0.083); 
this procedure results in 1512 measured bins. 

The quadratic estimator does produce a Fisher matrix as an estimate of error covariance 
of the measured data-points. In our case, this assumes that the flux fluctuations correlation 
within individual quasars dominate the covariance matrix. In practice we would like to have 
a more robust measurement of the error-bars for our final measurement; for this purpose we 
use various techniques that estimate the covariances internally. 

3.4 Distortion of the three-dimensional correlation function 

It has been known since the 1970s [51] that the naive estimates of 3D correlations from the 
Lyman-a forest result in distorted correlation function due to continuum fitting. 

The effect appears because continuum fitting inevitably entails inferring information 
about continuum shape from the actual forest data and the model will naturally try to min- 
imize the variance of residuals with respect to the predicted mean continuum. For example, 
an unbiased estimates of ^i^^ will have a small but finite mean due to presence of Fourier 
mode whose wavelength is much larger than then the length of the forest. Continuum fitting 
cannot distinguish between a true underlying offset and a change in quasar amplitude and 
will typically set the amplitude of such mode to zero. The size of this effect is surprisingly 
large and will result in a biased correlation function measurement. 

In this work, we would, in principle, not need to worry about distortion effects, because 
we are trying to find an isolated bump in the correlation function, while distortion is a 



-16- 



smooth function that could be absorbed with a sufficiently flexible broadband. Nevertheless, 
prevention is always preferable to cure. 

There are several ways to address this distortion. The simplest one, used in [34], is to 
assume that the effect of continuum fitting is to force a set of measured dp values to have 
zero mean in each quasar line of sight, by subtracting its mean dp dp — {Sp). One can 
propagate this assumption through a simple analytical model and obtain a formula for the 
correction. Resulting expression is approximate, but worked well for the signal-to-noise used 
in the above paper. In this work we can do better than that by removing the effect of the 
continuum fitting by down-weighting linear combinations of datapoints that we know are 
affected by the continuum fitting process, namely the constant offset and slope in per-quasar 
data-vectors. This is done by associating large variances with these particular modes as 
described in Appendix B. This technique is equivalent to discarding a point by associating 
a very large variance with it, only that we now associate a large variance with a particular 
linear combination of points. Of course, the contribution from a particular linear combination 
can only be down- weighted if one uses matrix weighting. The result is that, by construction, 
these marginalized modes cannot have an effect on the final result. 

In particular, we marginalize modes that appear as a constant offset and those that 
appear as a linear function of log A in the 6^ field. We add these terms to the covariance 
matrix both for the ID P(k) measurements as well as the 3D correlation function.. 

As a result of this marginalization, the estimator propagates the uncertainty associated 
with the unknown modes into the measurements of the 3D correlation function. As an 
unfortunate consequence of working in configuration space, this approach results in large, but 
nearly completely correlated errors in measurements. In Fourier space, this result corresponds 
to a few poorly measured low-A: modes. 

For plotting purposes, these modes can be self-consistently projected out of the mea- 
sured correlation function and theory predictions, resulting in a distorted correlation function 
measurement whose error- covariance is considerably more diagonal as we will do in the Sec- 
tion 6. However, the shape of the distortion is now known from the eigen-vectors of the 
projected modes. If one limits the analysis to inspection of values, this is, of course, 
unnecessary. 

3.5 Determining the covariance matrix of the ^ measurements 

There are two aspects of any analysis that can affect the results in a detrimental way: the bias 
of the estimation procedure and the uncertainty in the derivation of the covariance matrix 
of the errors. The former can be most precisely determined by examining the results on the 
synthetic data. As far as accuracy of the covariance matrix is concerned, fifteen realizations 
is not enough if one wants to consider events that are rarer than 1 in 15 (i.e., > 2a events). 
Therefore one must rely on the internal checks of the data, and we provide several such 
methods. 

3.5.1 Bootstrapped covariance matrix 

We start by creating a bootstrapped covariance matrix. This process recreates bootstrap 
"realizations" of the data by randomly selecting, with replacement, N plates from the original 
set of plates and combining them using their estimator provided weights. In each realizations 
some plates appear more than once and sometimes never. The covariance matrix of these 
realizations can be assumed to be a good approximation to the covariance of the estimation, if 
the plates are sufficiently independent and if we have sufficient number of them. Although the 
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plates are not strictly independent (the pairs of quasars spanning two plates are randomly 
assigned to one of the plates), the approximation that these are independent should be a 
very good one - as we will see later, not just plates, but even quasars can be considered 
independent. There is one subtlety, however. The variance associated with the distortion of 
the correlation function described in Section 3.4 will not be accurately represented in this 
bootstrapped matrix, because this variance will be artificially low (by virtue of these modes 
being calibrated out) and therefore the bootstrapped matrix cannot be used directly on the 
measured correlation function. 

The full bootstrapped covariance matrix from the dataset (or one realization of the 
synthetic data) is non-positive definite, because its size is 1512 x 1512 while the number of 
plates contributing to the bootstrapped matrix is only 817. In other words, the 817 vectors 
do not span the full 1512 dimensional space and therefore linear combinations of them cannot 
generate 1512 independent eigen- vectors of the covariance matrix. However, the bootstrapped 
covariance matrix where elements corresponding to different angular separations are assumed 
independent is positive-definite as we are now determining eighteen 84 x 84 matrices from 
eighteen sets of 817 vectors of size 84. We show that this approach is a good approximation 
in the next subsection. 

The bootstrap matrices used here were created from 150,000 bootstrap samples, which 
is well over what is required for convergence (covariance matrix converges in the sense of 
values changing at less than unity at about 30,000 samples). 

3.5.2 Comparing bootstrapped and estimator covariance matrices. 

The covariance matrix derived from our estimator has, by construction, structure that is 
diagonal in the angular separation. This can be understood as follows: our weighting matrix 
(discussed in Section 3.3) approximates quasars as being independent. This is the same as 
assuming that any given quasar pair is completely independent from any other quasar pair. 
Since a quasar pair contributes to the correlation function at only one angular separation 
bin (i.e. the angular separation between quasar lines of sight), it necessarily implies that 
different separation bins have uncorrelated errors. This assumption is not exactly true and 
we need to investigate the accuracy of this approximation. 

To this end we combined bootstrap matrices from 15 realizations of the mock datasets, 
i.e. from the full 15 x 817 = 12225 plates. These now hold enough information to create 
a fully positive-definite covariance matrix by bootstrap. We then inspect the correlation 
elements in this matrix that span different separations. These elements are small, typically 
less than 0.05. We plot the cross-correlation matrix in the angular separation direction for 
a fixed choice of 5 log A and separation bins in Figure 8. While the eye can discern some 
coherent structures, they do not dominate the matrix. 

To make this statement more precise we turn to fitting the bias parameters. We fit 
for parameters with their redshift evolution with the full bootstrapped covariance matrix 
and the one in which we have set the non-diagonal elements at different separations to zero. 
The fitted parameters were unchanged, and the best fit changes by less than 10 units 
for combined 15 realizations, implying a change of less than unity in per realization. 
We therefore conclude that for the present analysis, it is a good approximation to assume 
correlation function measurement errors to be uncorrelated across angular separations. 

Next we focus on the correlations inside a single angular separation bin. The quadratic 
estimator covariance matrix has two large eigenvalues corresponding to marginalization over 
amplitude and spectral slope that were performed as described in Section 3.3 and Appendix 
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Figure 8. Cross-correlation coefficients for the bootstrapped matrix in the separation direction, where 
other two parameters are held at fixed position in the middle of our fitting range ((5 log A = 0.023, 
z — 2.5). The diagonal elements, which are unity by definition, have been set to zero to emphasize 
the off-diagonal structure. 

B. The bootstrapped covariance matrix of course does not know about this marginalisation, 
so we manually remove these from the quadratic estimator covariance matrix. This is done 
by diagonalizing the matrix and then reconstructing it with the two largest eigenvalues set 
to zero. The resulting correlation matrices are displayed in Figure 9. This result shows that 
the basic structure of the correlation matrix is correct. However, it is difficult to judge these 
matrices in detail from the plots and therefore we turn to another test. 

3.5.3 The SK test 

We measure the correlation function in terms of plates. The quadratic estimator provides 
a covariance matrix estimate for each individual plate. These are then combined into the 
complete estimate using 

C^otito.^Y.'^-Hi (3.5) 

Ctol = J2cr' (3.6) 

i 
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Figure 9. Bootstrapped correlation matrix at a fixed separation (120^ but plots for other separations 
are similar). The left panel shows the bootstrapped covariance matrix, while the right panel shows 
the quadratic estimator matrix. Three 3x3 "tiles" in each plot correspond to the three redshift bins, 
while structure inside each "tile" correspond to the correlations between J log A bins. See text for 
further discussion. 

where the subscript "tot" denotes the total estimate and the subscript i the contribution 
arising from plate i (note that i does not denote the vector element). One could make an 
approximation that within errors of a single plate, the survey mean is a good approximation 
for the true mean. However, the resulting would be lower than expected, because any 
given plate has contributed to the total mean. However, it can be shown that the following 
pseudo [36] 

Axf = - ^totf {Ci - Ctot)-' i^i - Ctot) (3.7) 

has the expected mean for the given number of degrees of freedom. Additionally, if the 
distribution of the ^'s is Gaussian then these quantities are also distributed. 

We have performed this test on three sets of data: i) on the reduction where the con- 
tinuum and other nuisance quantities were assumed to be known perfectly, ii) on the fully 
analyzed synthetic data and iii) on the data. Results of these tests are presented in Figure 
10. 

From these plots we can make the following conclusions. First, even when all nuisance 
quantities are known, the estimator covariance matrix is not perfect as it underestimates the 
variance by about 10%. This can be due to several reasons, such as the ID power spectrum 
not being well estimated, but most likely arises from the fact that the field itself is non- 
Gaussian. However, the mis-estimate does not seem to correlate with separation, redshift or 
the magnitude of the eigenvalue contributing to it. When we move our attention to a more 
realistic analysis in which we fit for the continuum, etc, we find that the estimators make 
an overall underestimate of the variance by about 35%, but also that this underestimate 
is primarily produced by the large eigenvalues of {Ci — Ctot)- In other words, the noisier 
eigen- modes of the covariance matrix are more underestimated than the less noisy ones. We 
find no structure if we slice the data in other parameters (separation, redshift, etc.) 

This finding suggests an alternative way of correcting the covariance matrix. We fit the 
ratio of measured to expected scatter in Figure 10 as a function of eigenvalue rank with a 
smooth curve and apply it to each covariance matrix separately. The so-constructed matrix 
will have desired properties in terms of the SK test by construction. 



-20- 



200 
^ 150 

; 100 





200 
^ 150 
S 100 



0.8 0.9 1.0 1.1 1.2 1.3 1.4 
Plate <x^>ldoi relative to combined 



1.3 
1.2 


























1.1 


















1.0 










0.9 








■ 


0.8 







200 


400 


600 800 



r i.o; 
;^o.8: 

\ 0.6; 

\ OA 
0.0: 



200 400 600 800 1000 1200 1400 
Eigenmode Rank 



1.0 1.2 1.4 1.6 

Plate <x^>ldof relative to combined 



3.5 
^ 3.0 
JC2.5 
v^2.0 
I 1.5 

a 0.5 
0.0 







































200 400 600 
Plate Index 


800 
























































200 400 600 800 1000 1200 1400 
Eigenmode Rank 



0.9 1.0 1.1 1.2 1.3 1.4 1.5 
Plate <x^>ldof relative to combined 



A- 1.0 































)-, , , 








200 400 600 800 
Plate Index 








































































200 400 600 800 1000 1200 1400 
Eigenmode Rank 



Figure 10. Result of SK test done on the pure flux field in one realization of synthetic data (left 
column), on realization of the fully reduced synthetic data (middle column) and the real data (right 
column). For each dataset we plot the distribution of the per d.o.f. values (top row), the P^r 
unit d.o.f. as a function of plate number (middle row) and the relative contributions to this quantity 
coming from eigenvalues of Ci — Ctot sorted by eigenvalue rank and averaged over all plates (bottom 
row). See text for discussion. 



3.5.4 Choice of covariance matrix 

To recap this section, we have found that the estimator-provided covariance matrix is a 
decent approximation, but not sufficiently accurate to provide reliable error-bars. Hence we 
used three different methods for estimating the errors on our BAO parameters: 

• Method 1. Perform the eigen- value decomposition of the estimator matrix into eigen- 
vectors and eigenvalues cr|, and for each eigen- vector we determined the variance 
predicted by bootstrapped covariance matrix: 

^^>s = vfCbsV^ (3.8) 

We set the variance to the larger of the two: cr| max(a|, cr|^g) and reconstruct the 
covariance matrix. 

• Method 2. Use the eigenvalue-ranked fixed individual matrices as described in Section 
3.5.3 

• Method 3. Use the bootstrap technique on the final BAO parameters. Generate boot- 
strap samples as described in Section 3.5.1, then proceed to fit the BAO parameters 
with it. Use the final distribution of these parameters as an estimate of the uncertainty 
in those parameters. 
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As we will see, these estimates provide compatible results in the vicinity of high- 
likelihood regions. They do, however, differ in the tails of distributions; we will return 
to this point later. 

It might be argued that one should take the version of Method 1 in which the roles of 
bootstrap and estimator covariance matrices are reverse, i.e. bootstrap matrix is expanded 
eigenvalue-by-eigenvalue and small eigenvalues corrected by their variances in the estimator 
matrix. The problem with this approach is that in general eigen- vectors of the bootstrapped 
matrix are not orthogonal to the 36 large-eigenvalue eigen- vectors in the estimator matrix. 
We attempted various ways to address this issue, but none produced a satisfactory matrix (i.e. 
one that would pass tests described in Section 4). We are therefore omitting this approach 
in the paper. 



4 Fitting Baryonic Acoustic Oscillations parameters 

The process of fitting the BAO oscillation parameters is described in detail in the companion 
paper [36]. In short, we model the data as 

Cobserved(AlogA,5,2:) = Ccosmo (^|| , ^± , «|| , (1 + S^(r, ^, z)) + ^^(r,^,^). (4.1) 

That is, using the standard cosmological model, we convert the observed coordinates (A log A, 
5, z; difference in logarithm of observed wavelength, separation and redshift) into cosmological 
coordinates (r, fi, z) (radial distance in correlation function, angle with respect to the radial 
direction and redshift). Functions Br^ and Ba are multiplicative and additive broadband 
model discussed below. Our model for ^cosmo is given by 

Ccosmo = Cno peak(^||,^±) + ^peak ' Cpeak(^|| «|| , ^±Q^±) , (4.2) 

where we have decomposed the linear theory correlation function into a smooth component 
and the peak, but as opposed to most other BAO work, we only dilate the peak part of the 
correlation function. As discussed in [36], this is the most conservative way of ensuring no 
information is received from the broadband shape of the correlation function. The cosmolog- 
ical model has an adjustable "peak amplitude" parameter apeak- This amplitude is usually 
set to unity (recovering the full model), unless we explicitly state that we are fitting for it. 

The two isotropic dilation parameters measure the ratio of the BAO standard ruler size, 
namely sound speed r^, to the appropriate distance scale 

with respect to our fiducial model, which is flat ACDM with f^^/i^ = 0.0227, 9.m = 0.27 and 
h — 0.7. In the above equations H{z) is the Hubble parameter and Da{z) is the comoving 
angular diameter distance. 

For most of the test, we perform fits with the isotropic dilation factor a^^o which corre- 
sponds to the fits where we set a\\ — [— o^iso)- These isotropic fits fit a single parameter 
and hence have higher signal-to-noise and are similar to fits to the monopole of the correla- 
tion function often performed in the case of galaxy BAO. The reason why our measurement 
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Name 


Type 


i 


e 


k 


# free parameters 


BBl 


Additive 


0...2 


0,2,4 





9 


BB2 


Additive 


-2. ..0 


0,2,4 





9 


BBS 


Additive 


0...2 


0,2 


0,1 


12 


BB4 


Multiplicative 


0...2 


0,2,4 





8 


BB5 


Multiplicative 


0...2 


0,2 


0,1 


10 


BB6 


Additive + Multiplicative 


0,1 


0,2,4 





11 



Table 1. Broadband models used in this work. This table shows span of power indices i/ and k used 
for different broadband model defined in Equation 4.5. When counting free parameters, one must 
take into account that some are perfectly degenerate with standard bias parameters (i.e. monopole in 
multiplicative broadband is equivalent to bias change). See text for discussion. 

cannot be usefully compressed to Dy = [((1 + z)Da)'^{cz/H{z))]^/^ quantity lies in the fact 
that the redshift-space distortion parameter /3 is considerably larger in the case of Lyman-o; 
forest, considerably increasing the signal-to- noise ratio in the radial direction. One cannot 
therefore assume that the relative amount of information coming from Da{z) and H~^{z) 
will simply follow the geometrical factors of 2/3 and 1/3 respectively. 

When fitting the synthetic data we use a linear model, because this is what the synthetic 
data assumed. However, when fitting the real data, we use the model that has been smoothed 
to model the weakly non-linear dark matter power spectrum at the redshift of interest. 
Strictly speaking, this approach is not necessary as it negligibly affects our ^ind only 
marginally increases our error-bars, as we discuss further in section 5.2.1. We assume that 
the interpolation in the estimate of the correlation function takes care of the finite bin-size in 
the correlation function measurement. The linear cosmological model at a given redshift has 
two parameters: the bias h and the redshift-space distortion parameter /3. Following [34], we 
use the independent parameters /3 and (1 + 

We also use a parameter that governs the redshift evolution of the clustering amplitude, 
75 = (ilog(6^^^)/dlog(l + z), where g{z) is the growth factor. This means that we never 
explicitly use the growth-factor: we calculate the linear prediction at the reference redshift 
and then scale the linear power spectrum (together with non-linear corrections) using 75. 
Because the inclusion of broadband parameters significantly degrades our ability to measure 
these parameters we institute a broad Gaussian prior on /3 = 1.4 zb 0.4 and 75 = 3.8 =b 1. 

We additionally fit for a multiplicative broadband (S^) and an additive one {Ba)- 
These two broadband functions are smooth and are designed to remove any non-peak BAO 
information as described [36] . They are defined to be of the form 

where is the Legendre polynomial of the ^-th order and where we fit around 10 
parameters simultaneously. In Table 1 we tabulate which combinations of i,^. A: terms we 
float in the 6 broadband models considered in this paper. We have also investigated less 
physically motivated broadband models such as those with odd powers of ji and stronger 
angular dependencies. All models gave consistent results. For more information see the 
companion paper [36] . 

Our fit uses data in the range 50/i~^Mpc — 190/i~^Mpc. This choice contains sufficiently 
large buffers at both ends of the fitting range to be able to determine the broadband param- 
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Broadband 


BBl 


BB2 


BB3 


BB4 


BB5 


BB6 


Best fit Mel 
Best fit x^ Me2 


904 ± 42 
908 ± 40.6 


904 ± 42.3 
908 ± 41.2 


898 ± 39.5 

899 ± 36.9 


904 ± 42.4 
908 ± 40.8 


901 ± 41.3 
903 ± 38.3 


899 ± 42.5 
902 ± 41.2 


100 X (ai3o - 1) 
Comb. r. Mel 
Comb. r. Me2 


-0.795 ± 0.749 
-0.677 ± 0.727 


-0.544 ± 0.761 
-0.303 ± 0.735 


-0.641 ± 0.748 
-0.514 ± 0.726 


-0.818 ± 0.705 
-0.551 ± 0.675 


-0.862 ± 0.675 
-0.518 ± 0.648 


-1.06 ± 0.728 
-0.905 ± 0.684 



Table 2. Mean and variance (printed as error-bar) of best fit fo^" different choices of covariance 
matrix (labels Mel and Me2 refer to Method 1 and Method 2 covariance matrix) and broadband 
model (BBl. . . BB6) for our 15 synthetic datasets. The effective number of degrees of freedom is 894 
for BBl. We also show the result of fitting for the BAO position when all 15 Ax^ contributions are 
summed. 

eters while preventing the sharp upturn in the correlation function at smaller separations 
from affecting our fit. We also always cut the data at A log A < 0.003 as these data are the 
most affected by the sky subtraction residuals [34] . 

For each model, we minimize varying all parameters except the parameter of interest, 
e.g., Q^iso oi" and then gridding the minimum over the parameters of interest. This 

procedure is equivalent to marginalizing over all the other parameters in the limit of those 
marginal likelihoods being Gaussian. This approach is a good approximation when the 
broadband parameters are not completely degenerate. We have found that combining all 
broadband models into general broadband models (with ^100 parameters) results in the 
constraints that are very degenerate (in the sense that the data cannot distinguish between 
terms in the broadband model). In such case, our explicit Monte Carlo tests have shown that 
the prior volume completely dominates the posterior, or equivalently, the Ax^ distribution 
in minimization becomes very non distributed (in the sense that incorrect solutions can 
be found that are up to Ax^ > 25 away from the true solution). 

We calculate uncertainty on the parameters using a Bayesian approach: we convert Ax^ 
values into probabilities and then measure the median and relevant points in the cumulative 
probability distribution. This is an important point as the distribution is non-Gaussian with 
the confidence limits being asymmetric at a 10% level for our data. 

4.1 Tests on synthetic data. 

We start by discussing the goodness of fit for our basic fit. These results are shown in Table 
2. This table shows that the 15 synthetic mocks on average produce a good x^ that is 
independent of the broadband model used. We show results for both Method 1 and Method 
2 covariance matrices. Table 2 also presents the results of combining fits to the isotropic 
BAO from all mock realizations; results are consistent with no bias in the method and are 
consistent between the two methods for estimating the covariance matrix. There is weak 
evidence that Method 2 might be performing slightly better than Method 1. 

Figure 11 displays the x^ contours for different broadbands for our synthetic data for the 
isotropic fit for the Method 1 covariance matrix. We also show the results of the combined 
fit, which was calculated by summing Ax^ contributions from each individual measurement. 
The curves are parabolic in the vicinity of the best fit position, but far from this point, they 
yield a constant for the additive broadband. In those cases, the inferred bias is always close 
to zero: the model decides that the penalty for having a peak at the wrong position is too 
large and instead replaces most of the correlations with the broadband model which then 
makes the fit independent of a[so. This is not the case for broadbands with a multiplicative 
component, and in this case we indeed see structure outside the most likely region. We also 
find that for some "lucky" realizations, all broadbands agree on the position of the peak and 
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have the same error-bars, while for some the differences can be staggering. In particular, 
multiplicative broadbands often latch onto noise features. 

The ability of the broadband model to completely replace the true cosmological model 
when fitting the peak at an undesirable position is unphysical. In other words, one should be 
able to use the information that there is no evidence for a peak at certain dilations, not just 
the evidence that there is a peak at others. We can cure this problem by instituting a weak 
prior on (1 + /3)b = 0.336 =b 0.12. This prior essentially states that at two sigma, at least one 
third of the total signal must be cosmological in origin. Note that while the central value of 
this prior is the same as the value measured in [34], the error is ten times larger, i.e. this is 
a really weak prior compared to what measurements say about this quantity, curves for 
those fits are plotted as dotted lines in the Figure 11. The bottom line is that outside the 
favored region, the Ax^ difference increases further, but in the region of high-likelihood it is 
not affected, so this prior will not be used for the final results. 

Results for Method 2 are quite similar to those of Method 1. In Figure 11 we also plot 
comparisons between methods for two ways of estimating covariance method (dashed vs solid 
lines). We see that Method 2 in general gives higher significance, but not uniformly so. 

Next we investigate the probability that ai^o > 1.0 on each individual mock. If our 
estimates are unbiased and errors correctly estimated, this probability should be uniformly 
distributed between and 1. We plot the cumulative distribution of this quantity in Figure 
12. 

This figure deserves some attention. From the figure it appears that there might be a 
mild tension with the expected cumulative distribution. The minimum and maximum values 
of p{aiso > 1) are around 2-3% and 93-99%, which are consistent with no bias. Plotting the 
distribution withp(aiso > 0.995) or p(aiso > 1.005) produces outliers in the 1 percentile range. 
Our synthetic datasets share the same continuum realizations and are thus not completely 
independent, with continuum errors not canceling out when one sums realizations. This 
results in any one realization being consistent with the correct solution (i.e., no outliers), but 
they are somewhat correlated leading to skewed plots in Figure 12. 

In the absence of a larger number of and more independent synthetic datasets, we 
associate a 0.5% systematic uncertainty with our BAO peak position fitting. 

Next we move to the anisotropic fit and relax the assumption of a± = q;||. This is 
not a purely statistical exercise as the two quantities probe different underlying physical 
parameters: a± is a measure of the comoving angular diameter distance to the effective 
redshift, while q:|| is the measure of the local expansion rate at the redshift of interest. If one 
wants to use our results to estimate cosmological parameters, it is important to distinguish 
between the two parameters. 

We begin by plotting the equivalent of Figure 11 for the anisotropic fit. This is presented 
in Figure 13. We show plots for one broadband (BB3) and the method 2 covariance matrix. 
Plots for the method 1 covariance matrix look very similar with larger error-bars. Plots for 
other broadband models again are similar, although the contours do move for realizations 
where BAO is only weakly detected. For the strong detections and for the combined fit, the 
contours are essentially the same. 

4.2 Determination of uncertainties 

It is clear that for realizations for which the detection of BAO is strong, the likelihood is 
Gaussian and the error-bars symmetric and Gaussian. In this work we calculate error-bars by 
applying a uniform prior of aiso between 0.8 and 1.2 and calculating the confidence limits from 



-25- 




Figure 11. Ax^ curves for 15 synthetic realizations of the fuh dataset and the results of combining 
those curves (bottom right panel). Different colors corresponds to different broadbands: red (BBl), 
green (BB2), blue (BB3), cyan (BB4), magenta (BBS), yellow (BB6). Dashed lines are the same for 
the Method 2 covariance matrix. Dotted lines are the same upon applying the bias prior. Note the 
different scale on the axes for the case of combined realizations. 

the appropriate points in the cumulative distribution function. We compare these results with 
errors derived from a Gaussian expansion of the likelihood around the maximum-likelihood 
point in Table 3. The 1-sigma error-bars are always adequately determined by the Gaussian 
approximation (although also systematically underestimated at a percent level). However, 
for 3-sigma error-bars, the Gaussian approximation can fail catastrophically (i.e. the 3-sigma 
error-bars are limited by the prior). This result is to be expected; when the significance 
of the BAO detection is low, the information becomes prior dominated at higher levels of 
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•D.O 0.2 0.4 0.6 0.8 1.0 ^-Q.O 0.2 0.4 0.6 0.8 1.0 

p(aiso) > 1-0) p(a,,o) > 1.0) 



Figure 12. The cumulative distribution of the fraction of reahzations with value of aiso above the 
fiducial value when fitting synthetic data. For correctly estimated errors, this distribution should be 
flat between zero and unity, producing a straight line for cumulative distribution, plotted as the thick 
black line. Colors correspond to different broadband models and are the same as in Figure 11. The 
left plot is for Method 1 covariance matrix while the right plot is for Method 2 covariance matrix. 
Faint gray lines are 100 realizations of 15 points drawn from uniform distribution. 



Set 


+Na 


-Na 


Na symmetrized 


Na Gaussian approx 


Re 7, N=l 


1.59 


1.58 


1.58 


1.53 


Re 7, N=2 


3.28 


3.22 


3.25 


3.06 


Re 7, N=3 


5.22 


5.07 


5.15 


4.59 


Re 8, N=l 


2.31 


2.27 


2.29 


2.19 


Re 8, N=2 


5.24 


4.94 


5.09 


4.39 


Re 8, N=3 


16 


15.7 


15.9 


6.58 


Re 12, N=l 


2.79 


2.59 


2.69 


2.45 


Re 12, N=2 


8.57 


5.88 


7.22 


4.9 


Re 12, N=3 


23.7 


14.2 


18.9 


7.35 


Re 13, N=l 


1.81 


1.76 


1.78 


1.74 


Re 13, N=2 


3.83 


3.59 


3.71 


3.48 


Re 13, N=3 


6.43 


5.79 


6.11 


5.22 



Table 3. Gaussianity of errors for select realizations for broadband model 2. Middle section corre- 
spond to Bayesian errors determined from cumulative PSF while the right column is the Gaussian 
expansion around the most likely point. 

confidence [52]. 

5 Results on Baryonic Acoustic Oscillation parameters 

We start with the goodness of fit for our data. Fixing BAO parameters to the fiducial cosmo- 
logical model and using no broadband, we can attempt to fit for the basic bias parameters. 
With the method 1 covariance matrix we find /3 = 1.7 it 0.9, 6(1 + /3) = —0.336 it 0.019 and 
75 = 2.7 =b 1.4 and — 896 with 941-3 degrees of freedom. With the method 2 covariance 
matrix we find ^ = 1.3 ± 0.5, 6(1 + /3) = 0.350 ± 0.016 and 75 = 2.3 ± 1.3 and = 958. The 
simplest model of Lyman-a forest flux fluctuations, namely that we are measuring a linearly 
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Figure 13. Ax^ contours for 15 synthetic realizations of the fuh dataset and the results of combining 
those curves (bottom right plot). Contours are plotted at probabilities that enclose 68%, 95% and 
99.7% probability. This plot is for broadband model 3 and the Method 2 covariance matrix. Color 
scale saturates at Ax^ = 25. 

biased tracer of the underlying dark-matter field with Kaiser-like redshift-space distortions, 
is thus a satisfactory fit with just three floating parameters and is consistent with [34] even 
when a completely different radial scale range is employed (10/i~^Mpc < r < 90/i~^Mpc in 
[34] and 50/i~^Mpc < r < 190/i~^Mpc in this work; lack of small scales is also the reason 
for increased error-bars in this work despite quadrupling of the number of quasars). Since 
we have not tested these measurements on synthetic data, we mention them here just to 
demonstrate that the values are not wildly incompatible with previous measurements. 

Table 4 is the equivalent of Table 2 for the data. All broadband models lead to an 
improvement in compared to the pure cosmology fit, beyond what is expected from extra 
degrees of freedom. This was not the case in the synthetic data. The values are low but 
not pathologically so compared to expectation. 

In Figure 14 we show the curves for the data. Even without the prior on bias (dotted 
lines), we have a solid 4-sigma detection even in the case of the most conservative BB3 prior. 
The same figure also displays how these constraints depend on the choice of covariance matrix 
(solid vs dashed lines). Perhaps more relevant are the probability plots, which are shown in 
Figures 15 and 16. Here we have converted the same x^ values shown before into probabilities 
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Broadband 


BBl 


BB2 


BBS 


BB4 


BB5 


BB6 


Mel 


855.3 


857.3 


830.7 


864.2 


849.6 


856.9 


Me2 


897.5 


900.2 


863.7 


908.4 


890.5 


900.1 



Table 4. Best fit values of for different ciioices of covariance matrix (labels Mel and Me2 refer 
to Method 1 and Method 2 covariance matrix) and broadband model (BBl. . . BB6) for data. Table 
2 but for data. See text for discussion. 



Broadband 


BBl 


BB2 


BBS 


BB4 


BB5 


BB6 


Mel 


-1.7 ± 2.06 


-1.62 ± 2.01 


-1.49 ± 2.02 


-1.04 ± 2.4 


-0.813 ± 2.49 


-1.19 ± 1.97 


Mel w 7^ 


-1.72 ± 2.04 


-1.62 ± 2 


-1.57 ± 2.01 


-1.1 ± 2.41 


-0.885 ± 2.47 


-1.24 ± 1.94 


Mel w prior 


-1.69 ± 2.1 


-1.63 ± 2.09 


-1.49 ± 2.04 


-1.04 ± 2.39 


-0.821 ± 2.47 


-1.17 ± 2.06 


Mel packed obs. 


X 


X 


X 


X 


X 


X 


Mel w Lin. th. 


-1.87 ± 1.86 


-1.88 ± 1.91 


-1.73 ± 1.78 


-1.53 ± 2.17 


-1.57 ± 2.07 


-1.78 ± 1.85 


Mel w alt. NL 


-1.6 ± 2.06 


-1.57 ± 2 


-1.45 ± 2.01 


-0.984 ± 2.38 


-0.716 ± 2.46 


-1.18 ± 1.95 


Me2 


-1.99 ± 1.9 


-1.75 ± 1.85 


-1.64 ± 1.89 


-1.55 ± 2.22 


-1.3 ± 2.25 


-1.4 ± 1.82 


Me2 w prior 


-1.98 ± 1.92 


-1.77 ± 1.91 


-1.64 ± 1.91 


-1.55 ± 2.21 


-1.3 ± 2.24 


-1.28 ± 1.9 


Bootstrap 


x 


X 


X 


X 


X 


X 


Bootstrap w prior 


X 


X 


X 


X 


X 


X 



Table 5. The best fit 100 x {aiso — 1) (that is percentage deviation from fiducial model) for various 
choices of broadbands and covariance matrices. 

and have also added the third method described in Section 3.5.4: bootstrapping on the final 
parameters. This test treats the bootstrap samples of the BAO parameters as if they were 
the MCMC samples, i.e., we inferred the probability by treating the number density of the 
samples as relative probability. Whether this is a statistically sound procedure is not clear 
as it mixes frequentist and Bayesian statistical approaches. Figure 15 demonstrates excellent 
agreement between the three methods for most broadband models, except for broadband 
model 5. The agreement likely stems from the fact that all statistical methods agree in the 
Gaussian limit. However, when we make the vertical axis logarithmic as we do in Figure 
16, we see that bootstrap samples have considerably more outliers than predicted by the 
Bayesian method. Perhaps the correct conclusion is that our error-bar estimates are sound 
only to some 3-a distance from the most likely point. 

These results are numerically condensed in Table 5. The bottom line is that the fit 
is quite stable with respect to the choice of covariance matrix and the broadband model. 
Results vary mostly by around ^ 1/4 and in a few cases by ^ 1/2 sigma. This result implies 
that choice of broadband model and covariance matrix are not crucial and that our dominant 
uncertainty is statistical rather than systematic. 

Next we turn to the anisotropic fit. The resulting contours for different broadband 
models are plotted in Figures 17 and 18. The data are consistent with our fiducial model 
and stable with respect to the broadband model used for up to two-sigma confidence limits. 
Beyond that limit, the choice of broadband significantly affects where the contours close. 
Therefore, given the current signal-to-noise ratio, it is impossible to give a broadband in- 
dependent statement for allowed range for confidence limits above 95%. We also see that 
constraints stemming from the method 2 covariance matrix are essentially the same in the 
shape, but tighter. 

5.1 Redshift dependence of BAO position and determination of Zeff 

We measure the correlation function in three redshift bins z — 2.0, 2.5 and 3.0. The esti- 
mator interpolates between these redshift bins and hence the redshift evolution information 
is not lost. In order to make cosmological inferences from our dataset, it is important to 
determine what is the effective redshift at which we measure the position of the BAO. 
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Figure 14. Variation of as a function of aiso for the data. Lines of different colors are for different 
broadband models, following the same color scheme as in Figure 11. Solid lines are for Method 1 
covariance matrix without applied prior. Dotted lines are upon applying the bias prior and dashed 
lines are for the Method 2 covariance matrix. 

To find the effect redshift, we add a parameter to our fit that governs the evolution of 
the BAO scale with redshift: — rflog«iso/rflog(l + In other words, we assume that 
the isotropic scale factor varies as 

<^Uz)=(-^Y\ (5.1) 

One can expect that measurements of o^iso and 7^, are correlated so that o^iso is best con- 
strained at the position where data are the most constraining. We illustrate this property in 
Figure 19, where we present the two-dimensional constraints on the o^iso-To; plane for three 
different choices of Zref- We see the expected behavior as the contours turn their direction 
as we move from a low reference redshift to a high reference redshift. However, contours 
are not completely Gaussian; this is expected, when the detection in any one redshift is not 
very strong. This feature also means that the precise effective redshift is a poorly defined 
quantity. Different broadbands produce tightest constraints on o^iso between z ^ 2.3 and 
z ^ 2.5. Therefore we set our Zeflf ^ 2.4. This should not matter for any realistic model that 
we might want to test using our data. If a cosmological model has a widely varying a across 
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Figure 15. Probabilities for isotropic parameter ai^o- The panels correspond to six broadband 
models. For each model we plot the Method 1 inferred probability with the solid black line, Method 
2 with the dashed red line and the bootstrap number density with the dotted blue line. 

small redshift range, then our approach of wide redshift bins fails. In such cases one would 
expect the peak to be smeared, but, as discussed in the following section, we do not detect 
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Figure 16. Same as Figure 15, but with a logarithmic vertical axis. 



any evidence for such effect. 

Finally, this figure also demonstrates that the parameter governing the redshift evolution 
of the BAO peak is consistent with zero, which is by itself a useful systematics check. 
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Figure 17. Contours for the anisotropic fit on the data for the Method 1 covariance matrix. Plots 
show different broadband models, starting with BBl (top left) to BB6 bottom right. Contours enclose 
68%, 95% and 99.7% probability. The color scale saturates at Ax^ = 25. 

5.2 Discussion of systematic effects 

The Lyman-o: forest correlation is a difficult field: systematic effects that affect the measure- 
ments of two-point functions are many and closely spaced. However, the BAO measurement 
is particularly robust. Once the existence of the peak is established in the data, there are 
few possible systematic effects that can affect its position. Perhaps a suitable analogy is 
measuring the redshift of a distant object: it is very difficult to spectro-photometrically cal- 
ibrate spectra even for a well designed modern telescope, but if line features are available 
in the target spectra, one can nevertheless measure the object's redshift with considerable 
precision. 

5.2.1 Non-linear effects on the correlation function. 

The dark matter field at the redshift relevant for our analysis {z ^ 2.4) is already weakly 
non-linear, although the effect is considerably smaller than that due to non-linearities (and 
scale dependent biasing) for galaxy surveys at considerably smaller redshifts. In particular, 
one expects ^ 6/i~^Mpc smoothing of the correlation function in the radial direction and 
^ 3/i~^Mpc smoothing in the transverse direction [36]. 
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Figure 18. Same as Figure 17 but for Method 2 covariance matrix. The contours at the same 
confidence level are tighter. 

We implemented this smoothing in an approximate manner as described in [36]. In 
general, the anisotropic smoothing generates multipoles beyond hexadecapole, that depends 
on both the smoothing parameters and the redshift-space distortion parameters. However, 
when the smoothing can be treated perturbatively, these multipoles can be ignored and the 
analysis restricted to ^ = 0, 2, 4 multipoles. 

Following [8] , we only smooth over the peak part of the correlation function, but we have 
shown explicitly that smoothing the entire linear correlation function affects negligibly. 
This is expected since we are using only distances r > 50/i~^Mpc. 

To further explore this effect, we have also tried a simple isotropic smoothing with kernel 
size of 3, 6 and 12 h~^Mpc. The error-bars and best-fit Table 6. Based on x^, we cannot 
distinguish the models. We also see that using a smoothed model increases the position error- 
bars at ^10% level and changes the best-fit peak position by a small fraction of standard 
deviation. Both effects are negligible, but we use the non-linear model for consistency. 

Finally, we note that using a model which is unphysically smooth (the 12/i~^Mpc 
isotropic smoothing) decreases the goodness of fit by ^ 7 units. Systematics in the data 
and processing would typically smooth the peak rather than sharpen it, suggesting that 
our interpolation scheme is sufficient and that the measured peak has the expected width. 
The sharpness of the measured peak also precludes fast variations of the peak position with 
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Figure 19. Two-dimensional constraints on the aiso-7a plane. We plot this only for broadbands BB2 
(top tow), BBS (middle row) and BB6 (bottom row) for reference redshifts z^ef = 2.2 (left column), 
Zref = 2.4 (middle column), Zref = 2.6 (right column). See text for further discussion. 

redshift, which would result in further smoothing of the BAO peak. 



Model 



linear 
non-linear 
Alt. non-linear 
3 h~^Mpc iso. smooth 
iso. smooth 
12 /i~^Mpc, iso. smooth 



Best Fit Ax^ 



831 
830.698 
830.73 
830.775 
831.29 
838.446 



_1 7Q + -L-«1 +3-78 +6.04 
-'-•''^-1.75 -3.52 -5.39 
_-! ^Q+2.05 +4.25 +6.77 

^•"^^-2 -4.02 -6.18 
_-i /|^+2.04 +4.21 +6.68 
-•-•^^-1.99 -4.01 -6.17 
7Q+1.89 +3.94 +6.27 
-'-•''^-1.82 -3.66 -5.59 

7+2.08 +4.3 +6.86 
-'-•'-2.03 -4.07 -6.24 
_1 nQ+2-95 +6.2 +9.99 
l.U^_2.88 -5.78 -8.86 



Table 6. Results of fitting smoothed and unsmoothed model for the Method 1 with broadband model 
3. See text for discussion. 



5.2.2 Data cuts 

We explicitly investigated a few possible systematic effects, by restricting the analysis to a 
subset of data and rerunning our correlation function analysis to see if the peak position is 
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affected more than one would anticipate given the amount of the data removed. In particular, 
we have performed the following tests: 

• Removing the data at the ultraviolet end of the spectrograph. We normally use the 
data redward of 3600A. The spectra at the ultraviolet end of the wavelength range 
are the least understood part of the data. They are difficult to calibrate spectro- 
photometrically due to nature of the calibration stars. This is exacerbated by the fact 
that quasar fibers were drilled in a plate position which, given atmospheric dispersion, 
maximizes the blue throughput relevant for the Lyman-a forest rather than the red 
throughput used for other objects and calibration stars. Hence for this test we remove 
the data between 3600Aand 3700A. 

• Spectro-photometric variability. This cut is another way of examining the effect of poor 
spectro-photometry at the blue-end of the spectrograph. For every quasar, we calculate 
two variability scores ("slope" and "normalization") as follows. For each individual 
exposure, the mean inverse variance weighted flux is calculated in two wavelength bands 
(4000-4500 Aand 5500-6000 A). The "slope" ("normalization") score is calculated using 
a weighted average of the deviations of the differences (sums) of the mean flux in the two 
bands for an individual exposure to the mean from the differences (sums) of the mean 
flux in the two bands for all the exposures corresponding to a single observation. We 
then remove 10% of the quasars that show the most exposure-to-exposure variability 
in the sum of "slope" and "normalization" variability scores. 

• A more conservative upper rest-frame cut. In our standard analysis, we use an aggres- 
sive upper rest-frame cut of 1210 A, going nearly up to the actual Lyman-a emission 
line. We check the results by restricting to a much more conservative 11 85 A 

• DLA cut. Normally, we use quasars that have an identified DLA system, but we exclude 
the data from 1.5 equivalent widths around the DLA. In this test we have removed all 
quasars that contained one or more DLA systems. 

• Balmer line cuts. We investigated the effect of cutting around the position of the 
observed frame Balmer lines. As discussed in the Section 3.1, we see features due to 
mis-calibration. If these mis-calibration features were correlated, these could produce 
correlated systematic signal at fixed radial distances. In order to check for the effect, 
we have reanalyzed our data with two separate cuts eliminating data in 4102 ± 15 Aand 
4341 ± 15A. 

The results of these tests are shown in Figure 20. None of the cuts affect the data 
significantly except potentially the 4102ABalmer line cut. This effect has been analyzed 
in [35] and we show explicitly that data are statistically consistent with the full sample in 
Appendix D. 

The same figure also shows the effect of adding a parameter describing the redshift 
evolution of redshift-space parameter /3, namely 7^3 = din /3 / d\n(l + z); this result also 
appears in the Table 5. This parameter has negligible effect, regardless of the broadband 
being employed. 

5.2.3 Other systematics 

There are additional potential sources of systematics which we argue are negligible in this 
work: 
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Figure 20. Results of the systematics test. The blue solid line is our default data reduction. We 
additionally plot results from the 3700Acut (blue dotted), variability cut (blue dashed), DLA cut 
(yellow dashed), Arest = 1185 A cut (yellow dotted) and Balmer line cuts (red dotted and dashed). 
Solid green line shows result upon addition of 7/3 parameter. All tests are done for the entire data for 
broadband model BB2 with method 1 covariance matrix. Other fits shows similar change. See text 
for discussion. 

• Metal contamination. Another important systematic effect might be metal contami- 
nation [53]. The strongest is Si III, which absorbs at 1206. 5A, producing a correlated 
absorption separated by 2271 km/s from the Lyman-a absorption by the same gas, 
corresponding to around 20/i~^Mpc of radial separation. The largest contribution is 
the Lyman-a -Si III contamination, which is suppressed by a factor of ^ 20 with re- 
spect to the Lyman-a - Lyman-o: correlation. In other words, the measured correlation 
function in the fiducial cosmology is the true cosmology plus a suppressed shadow with 
radial distances shifted by ^ 20/i~^Mpc. We thus expect the contribution from Si III 
to be negligible unless our amplitude accuracy begins to approach 0(^ 1/20) of the 
total peak height. This warrants further investigation along with the impact of weaker 
metal lines and we will do so in a future publication. 

• Binning artifacts. In principle, when one uses flat bins to estimate the correlation 
function, it is possible that the actual bin centers do not correspond to the true bin 
centers, thus skewing the measurement of the BAO position. For example, if the number 
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Model 



100 X (c^iso - 1) 



100 X {a\\ - 1) 



100 X {a± - 1) 



_-i 7+2.09 +4.36 +7.03 
-■-•'-2.04 -4.13 -6.42 
_1 QQ+1-91 +3.95 +6.23 
-•-•^^-1.88 -3.79 -5.83 
_-i ^20+2.04 +4.32 +7.42 
-'-•'^^-1.99 -4.11 -6.77 
_-, 7r+1.87 +3.9 +6.34 
-•-•'^-1.83 -3.73 -5.9 
_1 /IQ+2-05 +4.25 +6.77 

-L.^^^_2 -4.02 -6.18 
_1 A/1+1-91 +3-94 +6.2 

-'-•'^^-1.87 -3.77 -5.77 
_-, n/1+2-48 +5.46 +10.7 
-'-•'-'^-2.33 -4.72 -7.53 
_-i rr+2.28 +4.89 +8.55 
-L.OO_2,i6 -4.35 -6.77 
_n Q1Q+2.59 +5.93 +14.8 
U.010_2 38 _4 79 _7 55 

_-i 0+2.32 +5.01 +8.78 

-'-•'^-2.17 -4.35 -6.7 
_-, -iQ+1.95 +4.03 +6.48 
-'-•-'-^-1.98 -4.27 -17.4 

_-i /1+1.81 +3.72 +5.86 
-^•^-1.83 -3.83 -6.88 



_-i Q+3.58 +7.51 +11.9 
-•-•^-3.43 -6.88 -10.1 
_1 n5;^+3-18 +6.56 +10.3 

-'-•'-'^-3.11 -6.3 -9.63 
_-i q-i+3.51 +7.58 +12.3 

-'-•'^-'--3.3 -6.69 -10.2 
_n ^5;^Q+3-l +6.56 +10.6 

U.OO^_2 97 _6 03 -9.45 
-1 A^+3-46 +7.21 +11.4 
-^•'-*^-3.33 -6.67 -9.91 
_n <^98+3-14 +6.47 +10.1 

U.UZO_3 07 _6,2 -9.49 
_n 01/I+3.86 +8.29 +12.4 
U.Z1^_3 -7.23 -10.9 
n n'^Q9+3-4 +7.25 +11.5 
U.UOyz_3 25 -6.59 -10.2 
n 90^+4-38 +9.51 +12.5 
U.ZUD_3 9 _7 73 _;l1.5 
n QQ^+3-77 +8.11 +11.5 

_n A<^zL+3-73 +8.15 +12.2 
U.DD'^_3 42 _7 09 _ii.2 
PI -1 p:q+3.23 +7.01 +10.8 
U.loy_3 02 -6.21 -10.2 



_1 07+7-53 +17.4 +28.9 
-6.88 -14.2 -23.3 
_/L n^+6-47 +14.6 +27.1 

^•'-'^-5.95 -12.3 -19.9 
_9 9/1+7.44 +17.3 +29.7 

^•^^-7.06 -15.3 -25.2 
_^ 1+6.4 +14.4 +27.2 

^•-•--6.08 -12.9 -21.3 
_n Q^+7-39 +17.2 +28.7 

'-'•^^-6.71 -13.8 -22.4 
_o Q7+6.4 +14.5 +26.9 

'^•^'-5.92 -12.2 -19.8 
_o 90+10.8 +25.6 +32.7 

O.ZO_8 4 _iQ Q _24.6 
-f\ 1 ^+8-58 +22.4 +34.7 

'^•-'-^-6.84 -13.5 -20.5 
_o 1+9.4 +21.7 +31.9 
'^•-'--8.07 -16 -24.2 
-f\ 8A+7-35 +16.9 +30.5 

^•0^-QA9 -13 -19.8 
_9 9^+7.35 +17.7 +30.9 

Z.ZD_8 7g _24.5 _27.6 
_A Q/1+6.24 +14.6 +31.9 

^•^^-7.09 -21.3 -24.9 



BBl Mel 
BBl Me2 
BB2 Mel 
BB2 Me2 
BBS Mel 
BB3 Me2 
BB4 Mel 
BB4 Me2 
BB5 Mel 
BBS Me2 
BB6 Mel 
BB6 Me2 



Table 7. Marginalized 1,2,3-sigma confidence limits for different broadband models and covariance 
matrices. 



of pairs is a steeply falling/ increasing function of separation, then the actual mean 
separation of pairs contributing to a given bin of correlation function will be less/more 
than the nominal bin center. This is alleviated by linearly interpolating between any 
two bins in radial and redshift directions (instead of using top hat bins) . We have also 
tested our procedures on synthetic datasets with exactly the same geometry as the real 
dataset. 

• Wavelength and astrometric calibration. Radially, the BAO distance corresponds to a 
separation in wavelength. The relative wavelength calibration is better than 10~^ [32] 
and therefore well under the required precision. Transverse, the BAO distance corre- 
sponds to an angular separation of the order 1 deg. The angular separation between 
any two quasars is known to arc-second precision and in addition these errors are not 
coherent. This is another completely negligible effect. 

5.3 Consensus result 

We have presented a number of results which are reasonably consistent, but vary in precise 
values, because of different assumed covariance matrix and broadband model. However, it is 
useful to reduce this variety to a single consensus value. 

We show the marginalized confidence limits for all broadband models and both co- 
variance matrices in Table 7. These limits were derived from the cumulative distribution 
functions of marginalized likelihoods as described in the Section 4.2. The marginalized error- 
bars on the a± and ay are significantly larger than the errors on the isotropic fit. This 
degeneracy is due to large cross-correlation coefficient between the two parameters, which is 
estimated to be ^ —0.55 by fitting a Gaussian to the best-fit model. 

For the final results, we would like to quote a consensus result that captures the sta- 
tistical and systematic uncertainty. To this end we select the broadband 2 model with the 
less constraining (method 1) covariance matrix. Inspecting the table 7 we see that it has 
values and errors that are not at the edges of distributions. Figure 14 shows that it is the 
most conservative of the additive broadbands. We could use the most relaxing broadband 6, 
however, this model latches onto a noise feature at aiso ^ 0.8, which breaks its 3-sigma errors 
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assuming the noise feature is indeed not real (no other broadband can fit it and it becomes 
weaker with the more constraining method 2 covariance matrix). Therefore we absorb these 
uncertainties by specifying systematic errors that should contain them. Hence we quote our 
fit 

100 X (ciso - 1) = -1.6t|g t\l tit (stat.) ±1.0 (syst.) (5.2) 
for the isotropic fit and 

100 X (^11 - 1) = -1.3t|^ tl^j t^^l (stat.) ±2.0 (syst.) (5.3) 
100 X {a^ - 1) = -2.2t^;^ t\i (stat.) ±3.0 (syst.) (5.4) 

for the anisotropic fit. We are not quoting the marginalized 3-sigma error in the a^, as it 
is dominated by the top-hat prior on the Q:||-a^ plane. Note, however, that this is not true 
for the un-marginalized likelihood, as large parts of the a\\-a^ plane are excluded by our 
measurement, even for the less constraining models. 

The systematic error should be understood as having a top-hat probability distribution 
function, and thus while significantly affecting the 1-sigma confidence limits, its effect on the 
3-sigma confidence limits are considerably smaller. 

The fact that our quoted systematic errors are of the same order as one-sigma statistical 
errors is not coincidental. As the signal-to-noise ratio will improve with additional data, so 
will our ability to constrain the broadband model, and the differences between best-fit values 
for different broadband models will naturally decrease. 

5.4 Comparison to Busca et al. 

In [35] we presented an analysis of the BAO in the Lyman-a forest from a very similar dataset. 
The two analyses are complementary and the present analysis builds upon the previous work. 
The main methodological differences between the two can be summarized as follows: 

• The analysis in [35] uses more conservative cuts in their choice of data. They exclude 
all quasars in which one or more DLAs are present as indicated by the visual inspection, 
while we remove only 1.5 equivalent widths around the position of the DLA. We also use 
a slightly different DLA catalog from [43]. Moreover, we aggressively use as much of the 
forest as available (1036A-1210A rest-frame for non-BAL quasars), while [35] restrict 
the analysis to 1054A-1184A rest-frame. These different restrictions are responsible for 
most of the increase in the signal-to-noise in the present paper. 

• The continuum fitting procedure is different in the two analyses. 

• We use considerably finer binning in our correlation function. In particular we measure 
the correlation function in three redshift bins and for each redshift bin we measure it in 
18 transverse and 28 radial bins. The analysis in [35] uses 2D measurement as an inter- 
mediate measurement, but they eventually compress all data into 90 measurements of 
the monopole and quadrupole of the correlation function for a given fiducial cosmology. 
The finer binning has an advantage that it is in principle easier to catch systematic 
errors that are localized in particular bins of the measured correlation function; how- 
ever, the larger number of bins makes measuring and understanding of the covariance 
matrix considerably more difficult. 
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• The analysis in [35] uses diagonal weighting with a fixed assumption about the am- 
plitude evolution of the forest correlations with redshift, but allow for a completely 
unconstrained /3 parameter (we institute a weak prior on /3 to regularize fits) 

• The analysis in [35] makes stronger assumptions about the broadband and in particular 
do not consider a multiplicative broadband. 

• The analysis in [35] uses a wider fitting range 20 — 200/i~^Mpc vs 50 — 190/i~^Mpc 
employed in this paper. 

• The analysis in [35] does not take into account the broadening of the peak by non-linear 
evolution. We show that this has a negligible effect on the position and only small effect 
on the errors. 

In [35] the isotropic dilation factors has been determined to be 100 x (aiso — 1) = 1 it 3, 
to be compared with our result 100 x {aiso — 1) = — I.6I2 o -4 1 -6 8 (stat.) ±1.0 (syst.). This 
is a one sigma difference, from what might appear to be essentially the same dataset. We use 
looser cuts on data, but our detection significance has improved more than one would expect 
based on the increase in the number of data-points alone. The same has been observed, when 
more inclusive cuts have been applied to the [35] analysis (N. Busca, private communication), 
so this result seems to be a feature of the sample. Given the increase in the signal-to-noise 
ratio, the difference is not statistically significant. 

5.5 Consistency with SPT and BOSS galaxy BAO results 

Recently, the measurements of the acoustic oscillations in the Cosmic Microwave Background 
by the SPT experiment [54], combined with WMAP7 results [55] predicted a position for the 
BAO peak z — 0.57 in the simplest flat ACDM model, which is in tension with the 
measurements of the BAO in the CMASS galaxies from the BOSS experiment [8]. We briefly 
discuss this tension in the context of our measurement. 

We show results for the quantities of interest for several standard models in Table 8. In 
the cosmology of the CMASS fiducial model, the SPT result can be cast as a measurement 
ol a — 0.971 ± 0.020, while the "consensus value" (after reconstruction) for the CMASS 
measurement is a = 1.033 ± 0.017. The difference 8a = 0.062 is discrepant at 2.3-sigma 
assuming the errors can be added in quadrature. We will not discuss the details of this 
discrepancy (see e.g. [56]), but instead ask if Lyman-a forest BAO can add anything to this 
discussion. To this end we added two modification of the CMASS fiducial model to Table 8 
that have the same value of Vs/Dy as the CMASS best fit, either by changing the value of 9 

or CJdm. 

Inspecting Table 8 show that the SPT best and CMASS best models predict values of 
(ay, a^) that differ by less than our uncertainties. Differences in a\\ are of the order of 0.5%, 
while changes to are of the order of 2%, both small compared to our error-bars. When 
correctly taking into account the covariance between a\\ and a^, the size of the effect is larger, 
but BOSS Lyman-o: forest BAO is currently unable to distinguish between the models in a 
statistically significant manner. 

Therefore, Lyman-a forest BAO will not be able to distinguish between these models 
in the ACDM scenario. This is because the Kaiser redshift-space distortion parameter (3 is 
large for the Lyman-a forest, pushing the statistical weight of the BAO signal towards the 
measurements of Hubble parameter (where the effect is small) rather than angular diameter 
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distance (where the effect is large). This is an expected result [3]. Lyman-a forest BAO 
measurements do not improve the constraints to the minimal flat ACDM model. It does, 
however, help constrain models with non-zero curvature or dark energy that is dynamically 
significant at high redshift. 
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2.5544 


0.995 
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35738 
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Table 8. Cosmological models relevant for this work and the CMASS/SPT tension. Columns are 
predictions of various quantities for the selected models. After specifying uj^ and uj^m one needs to 
specify either the value of the Hubble parameter h or equivalently = Vg/Di^ss^ the ratio of the sound 
horizon, to the distance of the last scattering. The Lyman-a fid refers to the fiducial cosmology 
used in this work. SPT best is the best fit model coming from SPT+WMAP7. The CMASS fid is 
the fiducial model used in [8]. The a values in the z = 2A case are calculated with respect to the 
Lyman-a fid model (and can be compared to a constraints in this paper), while those for the z = 0.57 
model were calculated in the CMASS fid model (and can be compared to a constraints in the [8]). 
CMASS best 1 and 2 are changes to the CMASS fiducial model that reproduce the best-fit position 
of the measured BAO position by either changing 6 or c^dm 

6 Visualizing the peak 

Errors in the correlation function measurements are always significantly more correlated than 
the errors in a more natural space for two-point function such as Fourier space measurement. 
In our case, the marginalization over the unwanted modes makes this issue particularly 
acute - one obtains unbiased estimates at the price of large, but nearly perfectly correlated 
uncertainties. We also measure the correlation function in many points, each of which is 
rather noisy. Moreover, some of the BAO signal is in the monopole while some is in the 
quadrupole. Nevertheless, one would still wish to simply "see" the BAO peak and produce 
the peak plots in an objective manner. 

To this end we compress the measured correlation function as follows. We fit for a 
cosmological model with broadband. We then fix the parameters of this model, set the BAO 
amplitude to zero and fit for the binned deviations of the real space correlation function. 
The real space correlation function is piece-wise linearly interpolated between bins in the 
^(r)r^. Given the best-fit /3 value we propagate this fit to the monopole to higher multipoles. 
The result of this exercise is an optimally estimated real-space correlation function that has 
Kaiser-like redshift-space distortions. Since we are fitting around the model without the BAO 
peak, this fit has no way of knowing about the existence or the position of the BAO peak. 
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However, the resulting correlation function will still have large correlated errors. We can 
remove these, by diagonalizing the covariance matrix of our fitted binned correlation function 
and setting the amplitude of the largest eigen-mode to zero. At the same time, we must also 
project out the same eigen-vector from the theory. We have thus made a full circle; we end 
up with a distorted correlation function and a theory that accounts for this distortion. 

Results of this exercises are plotted in Figure 21 for data and for the mean of all 
realizations of synthetic data. Removing more than one eigenvalue produces essentially the 
same plots. 
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Figure 21. Correlation functions under the assumption of Kaiser-like redshift-space distortions for 
all fifteen realizations of synthetic data combined (left column) and the real data (right column). The 
top row are bins as inferred, the bottom row is with the largest eigenvalue removed. The dotted line 
is the model around which we measure the binned correlation function, the dashed-line is after taking 
into account the eigen- vectors projected out. The blue line is the theoretical peak added back to the 
dashed-line model to show agreement. The correlation functions are normalized to the underlying 
dark-matter correlation function with the bias factor inferred for our fiducial cosmological model 
divided out. 



7 Discussion &; Conclusions 

In this paper we have revisited the measurement of the Lyman-a forest BAO. We have 
detected the baryonic acoustic oscillations in the flux fluctuations in the Lyman-a forest of 
distant quasars which allows a measurement of the expansion rate at a redshift of 2: ^ 2.4. 
We have confirmed an earlier result with a significantly different method and data cuts. 



-42- 



Our results are consistent with but more stringent than [35]. By far the most important 
reason for this improvement is that we use more data by imposing looser cuts (which we 
think is valid when one tries to find an isolated feature such as BAO peak). Second, we use 
an optimal estimator, which should increase the available signal-to-noise ratio by 10-20% as 
indicated by tests on fields with known continuum. Finally, we use a much finer binning, 
resulting in less information loss when simultaneously fitting for the BAO and broadband 
parameters. 

Our results are more precise than the theoretical expectations for what the BOSS ex- 
periment with the present dataset is expected to achieve with the Lyman-a forest technique 
[25]. The covariance of and a\\ measurements have the same correlation coefficient as 
predicted by the Fisher matrix forecast —0.55), but our signal-to-noise is about 40% 
better. The Fisher matrix calculation assumed bias parameters from [19] and the forecasted 
expectations for the BOSS spectrograph signal-to-noise. In fact we have now achieved the 
statistical position accuracy that is consistent with what would be expected at the end of the 
survey. We saw with the synthetic data that realization-to-realization scatter is large so it is 
possible that we might just have been fortunate. On the other hand, it is also possible that 
the forecasts were pessimistic and that the signal is larger than expected. To distinguish 
between these possibilities, one needs to perform the analysis on a larger dataset and/or 
carefully compare the achieved spectrograph signal-to-noise and measured bias parameters 
with values that were used in the Fisher matrix analysis. 

We fitted different broadband models and found that these produce shifts in the best-fit 
values of BAO position of the order of the 1-sigma statistical errors. This effect is similar to 
the difference between the correlation function analysis and the power spectrum analysis in 
[8] - there is just one two-point function and hence different results reflect how a broadband 
model parametrized in Fourier or configuration space is simply a different model leading to 
different results. As the signal-to-noise increases, the likelihoods will become more Gaussian, 
which will make it easier to make sense of differing results from different broadband models. 

We attempted to be conservative when quoting our results: we used the least constrain- 
ing broadband model and the less constraining of our two covariance matrices. We also 
assigned differences between best-fit models from different broadbands to systematic error. 
Using these criteria, we measure 100 x {a^Qo — 1) = — 1.6^2 o ^4 1 ^6 8 (stat.) ±1.0 (syst.) with 
respect to the fiducial model (flat ACDM with = 0.27, h — 0.7) for the isotropic dila- 
tion factor and we find marginahzed constraints 100 x {a\\ — 1) = —1.3^33 tio 2 (stat.) 
±2.0 (syst.) and 100 x {a^ - 1) = -2.2^7;^ (stat.) ±3.0 (syst.) for the anisotropic fit. 

We do not attempt to present a precise number for the statistical significance of the 
peak detection. Assuming our most optimistic broadband model, our significance is well 
over 5-sigma. Believing that the main solution is the correct one, but choosing our most 
conservative broadband models reduce it to ^ 4-sigma. Taking seriously the secondary 
solution at a^^o ^ 0.8 found by the broadband model 6, the significance drops to ^ 3-sigma. 

To a large extent, the question of the BAO peak significance is a moot point. One is 
measuring the BAO position in order to learn about possible and impossible cosmological 
models. Plausible cosmological models will not deviate from concordance by more than a few 
percent and hence stay in the region where all methods give consistent results on the peak 
position. 

An important lesson learned in this work is that it is essential to have a method that is 
computationally feasible enough to analyze not just the real dataset, but many (hundreds) 
of realizations of the mock data to gain full understanding of the estimator and possible 
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systematics effects arising from the analyses. What we might have gained in statistical errors 
was lost in our ability to fully characterize the estimator using a large number of realizations. 
However, when going beyond measuring BAO and attempting to use the broadband power to 
constrain cosmology, an approach similar to the one employed in this work will be necessary 
to deal with the continuum-fitting induced distortion. In addition, a better analysis requires 
a more sophisticated continuum fitting and instrument model and will require a red-side 
systematics subtraction. In short, it is clear that many improvements can be made in the 
Lyman-o: forest data reduction methodology. 

We defer deriving cosmological constraints from the position of the BAO peak to a future 
publication. The current measurement of the BAO peak position is a striking confirmation 
of our basic understanding of the universe, as pointed out in [35]; and we expect it to be very 
competitive for models with early dark energy and curvature. 
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A Appendix: Optimal Quadratic Estimators 

Optimal quadratic estimators are a mature technology used extensively to estimate the two- 
point function of correlated data [47-49]. We will provide a general derivation in the next 
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subsection. This formulation is what we use to measure the ID power spectrum of quasars. 
However, for measuring the 3D power spectra, we use a shghtly modified procedure which 
halves the sizes of matrix operations in our case: an optimal cross-correlation estimator - 
this is discussed in subsection A. 2. 

A.l Optimal auto-correlation estimator 

Optimal quadratic estimators can be derived in several ways. The most common method is 
to write the Gaussian likelihood function, take its derivatives and then cast it as a Newton- 
Raphson approach to the minimum. There are subtleties to this derivation in that it produces 
the optimal unbiased estimator only if the second derivative is approximated as the Fisher 
matrix and not if one uses the exact second derivative in the Newton Raphson step. An 
alternative approach, which we follow in this paper, is to write a plausible ansatz for the 
estimator and then show in what limit it is optimal. 

Let d be a data vector of size whose correlation properties we wish to estimate. The 
expectation value of pairs of pixels is 



where C is the N x N covariance matrix, which we assume can be written as a sum over 
some component matrices C^^ = dC/dOi with coefficients 9i. We implicitly sum over repeated 
indices. If the covariance matrix is written in this manner, we can identify 9i as the estimate 
of the two-point function in bin i. In general, Oi can be values describing either the mea- 
surements of the correlation function or of the power-spectrum (depending in C^^). The J\f 
matrix is the noise matrix which we take to be provided and has no free parameters. 
We start by considering an estimator of the form 



where W is a weighting matrix which is at the moment a completely general symmetric 
matrix. 

The expectation value of Ei is given by 




(A.l) 



Ei = (Wd)^Ci(Wd), 



(A.2) 



{Ei) = Tr (CWC,iW) = GijOj + h 



(A.3) 



where 



hi = Tr (A/'WC,iW) 
Gi,- = Tr(C,iWCjW) 



(A.4) 
(A.5) 



One can therefore construct an unbiased estimator 



(A.6) 



where, by construction {^9ij — 9i. 

The (co) variance of this estimator is 





<(d'^WCfcWd-6fc)>((d'^WC;Wd-6;)> , (A.7) 
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Applying Wick's theorem to the average in the first term yields three terms. One of 
these terms decouples the first two brackets and thus cancels the second term and the other 
two are equivalent. One thus obtains 

Kij = 2Gr^^G-i^Tv (CWC,fcWCWC,iW) (A.8) 

Although these equations are complicated, they in fact reduce to the well-known cases 
for suitable choices of W, which we will discuss next. 

For simplicity, we will consider a toy model, where we try to measure the correlation 
function for a set of pixels. In this case C^i is unity for pairs of pixels corresponding to bin i 
and zero otherwise (in other words, we approximate correlation function as a binned function 
whose value is constant across each bin). Consider the case of a trivial estimator first. In 
that case W is an identity matrix and Ei becomes the sum over pixels pairs corresponding 
to correlation function in bin i. Gij becomes a diagonal matrix, whose elements are simply 
the number of pixel pairs corresponding to correlation function in bin i. Equation A. 6 thus 
reduces to 



pairs, i 



(A.9) 



In other words, we calculate the average of product of data pairs corresponding to the 
relevant bin - exactly as one would naively expect. What is the error of this estimator? 
Equation A.8 can be rewritten as 

Krj^^ E CacCw. (A.IO) 

-^^pairs,i-^^pairs,7' . . ■ 
pairs{a,b)&,{c,a)^j 

where sum is over all pairs (where a, b pair is counted as distinct from 6, a pair) of data-points 
with indices a and b that contribute to bin i and pairs of data-points with indices c and d 
that belong to bin j. The error estimation using this path is thus an operation, although 
the memory requirements are trivial. 

Optimal weighting can be obtained by setting W = C~^. In this case, the equation 
simplifies considerably. Both Gij and Kij can be expressed in term of Fisher matrix: 

F,, = Ixr (C,,C-iC,,C-i) = = Kr.\ (A.ll) 

Since the covariance matrix of estimator is the inverse of the Fisher matrix, the Crammer-Rao 
inequality stipulates that this must indeed be the optimal weighting. 
The estimator itself simplifies to 

^^ = lK^\E-bJ). (A.12) 

Often this inversion is unstable, especially if one is estimating bins that are highly-correlated. 
In that case, one can stabilize this inversion by instead citing the so-called filtered estimates 

m 

These estimates are now weighted with respect to the underlying valyes through the window 
function given by Pij)~^Pjk- For power spectrum measurement, these window functions 
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are bell shaped around the central value, but since the power spectrum varies smoothly with 
k and z, the resulting power spectrum estimate is essentially unbiased. 

Intuitively, the closer the weighting matrix is to the inverse covariance, the closer to 
the minimum variance our estimator becomes. The inverse variance weighting is just an 
approximation that the off-diagonal elements of the covariance matrix are negligible. 

In this work we use the same formalism for both ID and 3D two-point function mea- 
surement. For the ID power spectrum measurement, the matrices contain the response 
of the ID correlation function to a change in a single power spectrum band-power bin. We 
use the filtered estimate of Equation A. 13 and iterate until the measurement converges: at 
each iteration our weighting matrix uses the power spectrum estimated at the previous iter- 
ation. For the three-dimensional correlation function measurement, the matrices contain 
the response of the covariance matrix to a change in one correlation function bin (linearly 
interpolated in redshift) and we do not need to iterate: the weighting comes purely from the 
ID power spectrum measurement. 

A. 2 Approximating covariance matrix as block-diagonal 

The fully optimal quadratic estimator described above requires matrix operations on ma- 
trices of size N X which is impractical for the full Lyman-a forest dataset, where N = 
^quasars ^pixels ^ 10^. Howcvcr, while pixcls from a single quasar are very strongly correlated, 
they are only weakly correlated between quasars. Therefore it is sensible to approximate the 
covariance matrix as block diagonal, where each block corresponds to one quasar. In other 
words, we are assuming that quasars are independent as far as weighting is concerned. 

To make progress, let us collate the data vectors from the quasars into a single vector 
of size Nt = Ni + N2 (where Ni and N2 are the sizes of the datavectors from two quasars 
under consideration). The covariance matrix for this system can be written 



where we assume that the covariance matrix between pixels from a single quasar is known 
and we are only measuring the cross-correlations. matrices are of size A^i x N2 and so 
neither symmetric nor square. The derivatives with respect to the covariance matrix are then 
given by 



If the two lines of sight are so distant that the correlation between them can be neglected, 
then the inverse of the correlation matrix can be approximated as 



J:^C'i0^ C2 



(A.14) 




'0 c,- 



(A.15) 



.-1 ^ [cr^ 



(A.16) 




Plugging this result into optimal quadratic estimator, results in 



(A.17) 



where now we have 



(A.18) 
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with Fisher matrix 

Fii = = Tr (Cr^C,iC2-^C;5) . (A.19) 

There is no noise contribution, because we assume noise has no contribution that spans 
quasars (i.e. in the matrix of Eq. A. 14, the noise resides in the Ci and C2 matrices and is 
assumed to be diagonal in this work). 

These expressions are what one would naively expect for generalization of optimal 
quadratic auto-correlation equations. 

B Appendix: Marginalizing out unwanted modes 

To marginalize out templates, we use a standard technique from the Cosmic Microwave 
Background studies. Consider modifying the inverse of weighting matrix W by 

W'"^ = W"^ + Att^, (B.l) 

where A is a large number and t is the template we wish to marginalize over. It is clear that 
t is an eigen- vector of the matrix Att^ with an eigenvalue of Alt p. As A ^ oc it is true that 
t is also an eigen- vector of W'~^ and hence also an eigen- vector of with an eigenvalue of 

{A\t?)-'- 

Therefore, by rotating d into the eigen-space of the weighting matrix W, it is clear that 
the product W'(d + at) = W'd + a{A\t\^)-^ - Wd if A is sufficiently large. 

In other words, when dealing with a set of data-points, we can effectively discard a 
certain point by making its errors very large. Here we effectively do the same, but instead 
of discarding a single data-point, we discard a given linear mode by making its variance very 
large. This action is equivalent to a rotation into a space in which the offending mode would 
be a basis function, setting a large error on this one data-point that we want to discard and 
then rotating back. 

By performing the transformation in Equation (B.l), our estimator becomes insensitive 
to changes in the template, effectively marginalizing over it. The price paid, however, is the 
loss of signal-to-noise ratio in that linear mode. 

We wish to marginalize over the leading contributions from continuum fitting - the 
continuum amplitude and slope, namely t = const, and t = log(A/Ao). 

C Appendix: Coarse-graining 

When calculating the correlation function at large distances, it might prove useful to coarse 
grain neighboring pixels into larger pixels. To do this correctly, let us write the optimal linear 
estimator for the coarse-grained pixels. We will, for the moment, stop treating the field as 
random and write the likelihood for the coarse grained pixels as 

L = log£ = log(27r) - log \C\ SD)^C-Hd - SD)] , (C.l) 

where d are now the fine pixels in one given line of sight, are the values of coarse pixels and 

5 is the coarse-graining matrix which has dimension N x Nc, where N is the number of fine 
pixels and Nc is the (lower) number of coarser final pixels. For example, when coarse- graining 

6 pixels into 2, the matrix S would be of the form 
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'1 0" 

1 

01 
1 
.0 1. 

Taking the first and second derivatives of L with respect of Dj, we obtain 



(C.2) 



dL 
9D" 



= -[5'^C-Hd-SD)]^ 



[s^c-^s] 



(C.3) 
(C.4) 



The second derivative reveals that the covariance matrix of coarse pixels is given by the 
coarse-graining of the inverse of the fine-pixel matrix: 



Ccoarse " ^^^) — ^ 

Solving the first derivative for D, we find that 



Ccoarsel^ — S^C -^d 



(C.5) 
(C.6) 



In other words, the optimal coarse-graining is done by: 

• Calculating and C~^d 

• Coarse-graining both quantities to get C^^j.gg and C^^^-ggD 

• Inverting the coarse covariance matrix again and multiplying with the vector to get 
values of D. 

A trivial example is when assuming C = Sijaf. In that case, coarse graining reduces to 
the usual inverse-variance weighted mean. 



D Appendix: Statistical significance of the Balmer cut shift 

Figure 20 shows the effect of cutting the data in the vicinity of the Balmer 4102A line on the 
position of the BAO peak. How significant is this shift? 

In order to answer this question given our significantly non-Gaussian ? we subtract 
the cut-data Ax^ from the full Ax^. This difference is essentially the "pull" of the data that 
we cut to the total signal. The question is whether these are compatible. 

Assuming these two contributions to be independent, we can infer the probability that 
Q^iso of the two distributions are different. In order to do this, we infer the marginalized 
probability distribution for the difference between two aisQ. We plot results in Figure 22. The 
probability distribution for Aa has been calculated under the assumption that 0.8 < a < 1.2 
and does not go to zero away from the regions on high probability, because the cut data are 
consistent with any ai^o at a fixed Ax2 ^ 9. This plot shows that data are consistent with 
A^iso = at - 10% level. 
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Figure 22. Values of Ax^ for the cut data (red), complete data (blue) and difference (green) are 
plotted in the upper right corner. The implied probability distributions are plotted in the upper left 
corner. The inferred probability for Aaiso is plotted in the lower panel. See text for discussion. 
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